From 597b5ce028cee907e166afbb5aca7527ccb2d3e9 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 29 Jan 2024 17:22:27 +1100 Subject: [PATCH 1/5] migrate lectures --- lectures/_config.yml | 5 +- lectures/geom_series.md | 921 ---------------------------------------- lectures/lp_intro.md | 808 ----------------------------------- lectures/short_path.md | 487 --------------------- 4 files changed, 4 insertions(+), 2217 deletions(-) delete mode 100644 lectures/geom_series.md delete mode 100644 lectures/lp_intro.md delete mode 100644 lectures/short_path.md diff --git a/lectures/_config.yml b/lectures/_config.yml index 51848cf92..baf81c18f 100644 --- a/lectures/_config.yml +++ b/lectures/_config.yml @@ -102,8 +102,11 @@ sphinx: index_toc.md: intro.md # Remote Redirects redirects: - heavy_tails: https://intro.quantecon.org/heavy_tails.html + heavy_tails: https://intro.quantecon.org/heavy_tails.html ar1_processes: https://intro.quantecon.org/ar1_processes.html + geom_series: https://intro.quantecon.org/geom_series.html + lp_intro: https://intro.quantecon.org/lp_intro.html + short_path: https://intro.quantecon.org/short_path.html tojupyter_static_file_path: ["source/_static", "_static"] tojupyter_target_html: true tojupyter_urlpath: "https://python.quantecon.org/" diff --git a/lectures/geom_series.md b/lectures/geom_series.md deleted file mode 100644 index 1ffc0a706..000000000 --- a/lectures/geom_series.md +++ /dev/null @@ -1,921 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -(geom_series)= -```{raw} html -
- - QuantEcon - -
-``` - -```{index} single: python -``` - -# Geometric Series for Elementary Economics - -```{contents} Contents -:depth: 2 -``` - -## Overview - -The lecture describes important ideas in economics that use the mathematics of geometric series. - -Among these are - -- the Keynesian **multiplier** -- the money **multiplier** that prevails in fractional reserve banking - systems -- interest rates and present values of streams of payouts from assets - -(As we shall see below, the term **multiplier** comes down to meaning **sum of a convergent geometric series**) - -These and other applications prove the truth of the wise crack that - -```{epigraph} -"in economics, a little knowledge of geometric series goes a long way " -``` - -Below we'll use the following imports: - -```{code-cell} ipython -import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size -import numpy as np -import sympy as sym -from sympy import init_printing -from matplotlib import cm -from mpl_toolkits.mplot3d import Axes3D -``` - -## Key Formulas - -To start, let $c$ be a real number that lies strictly between -$-1$ and $1$. - -- We often write this as $c \in (-1,1)$. -- Here $(-1,1)$ denotes the collection of all real numbers that - are strictly less than $1$ and strictly greater than $-1$. -- The symbol $\in$ means *in* or *belongs to the set after the symbol*. - -We want to evaluate geometric series of two types -- infinite and finite. - -### Infinite Geometric Series - -The first type of geometric that interests us is the infinite series - -$$ -1 + c + c^2 + c^3 + \cdots -$$ - -Where $\cdots$ means that the series continues without end. - -The key formula is - -```{math} -:label: infinite - -1 + c + c^2 + c^3 + \cdots = \frac{1}{1 -c } -``` - -To prove key formula {eq}`infinite`, multiply both sides by $(1-c)$ and verify -that if $c \in (-1,1)$, then the outcome is the -equation $1 = 1$. - -### Finite Geometric Series - -The second series that interests us is the finite geometric series - -$$ -1 + c + c^2 + c^3 + \cdots + c^T -$$ - -where $T$ is a positive integer. - -The key formula here is - -$$ -1 + c + c^2 + c^3 + \cdots + c^T = \frac{1 - c^{T+1}}{1-c} -$$ - -**Remark:** The above formula works for any value of the scalar -$c$. We don't have to restrict $c$ to be in the -set $(-1,1)$. - -We now move on to describe some famous economic applications of -geometric series. - -## Example: The Money Multiplier in Fractional Reserve Banking - -In a fractional reserve banking system, banks hold only a fraction -$r \in (0,1)$ of cash behind each **deposit receipt** that they -issue - -* In recent times - - cash consists of pieces of paper issued by the government and - called dollars or pounds or $\ldots$ - - a *deposit* is a balance in a checking or savings account that - entitles the owner to ask the bank for immediate payment in cash -* When the UK and France and the US were on either a gold or silver - standard (before 1914, for example) - - cash was a gold or silver coin - - a *deposit receipt* was a *bank note* that the bank promised to - convert into gold or silver on demand; (sometimes it was also a - checking or savings account balance) - -Economists and financiers often define the **supply of money** as an -economy-wide sum of **cash** plus **deposits**. - -In a **fractional reserve banking system** (one in which the reserve -ratio $r$ satisfies $0 < r < 1$), **banks create money** by issuing deposits *backed* by fractional reserves plus loans that they make to their customers. - -A geometric series is a key tool for understanding how banks create -money (i.e., deposits) in a fractional reserve system. - -The geometric series formula {eq}`infinite` is at the heart of the classic model of the money creation process -- one that leads us to the celebrated -**money multiplier**. - -### A Simple Model - -There is a set of banks named $i = 0, 1, 2, \ldots$. - -Bank $i$'s loans $L_i$, deposits $D_i$, and -reserves $R_i$ must satisfy the balance sheet equation (because -**balance sheets balance**): - -```{math} -:label: balance - -L_i + R_i = D_i -``` - -The left side of the above equation is the sum of the bank's **assets**, -namely, the loans $L_i$ it has outstanding plus its reserves of -cash $R_i$. - -The right side records bank $i$'s liabilities, -namely, the deposits $D_i$ held by its depositors; these are -IOU's from the bank to its depositors in the form of either checking -accounts or savings accounts (or before 1914, bank notes issued by a -bank stating promises to redeem note for gold or silver on demand). - -Each bank $i$ sets its reserves to satisfy the equation - -```{math} -:label: reserves - -R_i = r D_i -``` - -where $r \in (0,1)$ is its **reserve-deposit ratio** or **reserve -ratio** for short - -- the reserve ratio is either set by a government or chosen by banks - for precautionary reasons - -Next we add a theory stating that bank $i+1$'s deposits depend -entirely on loans made by bank $i$, namely - -```{math} -:label: deposits - -D_{i+1} = L_i -``` - -Thus, we can think of the banks as being arranged along a line with -loans from bank $i$ being immediately deposited in $i+1$ - -- in this way, the debtors to bank $i$ become creditors of - bank $i+1$ - -Finally, we add an *initial condition* about an exogenous level of bank -$0$'s deposits - -$$ -D_0 \ \text{ is given exogenously} -$$ - -We can think of $D_0$ as being the amount of cash that a first -depositor put into the first bank in the system, bank number $i=0$. - -Now we do a little algebra. - -Combining equations {eq}`balance` and {eq}`reserves` tells us that - -```{math} -:label: fraction - -L_i = (1-r) D_i -``` - -This states that bank $i$ loans a fraction $(1-r)$ of its -deposits and keeps a fraction $r$ as cash reserves. - -Combining equation {eq}`fraction` with equation {eq}`deposits` tells us that - -$$ -D_{i+1} = (1-r) D_i \ \text{ for } i \geq 0 -$$ - -which implies that - -```{math} -:label: geomseries - -D_i = (1 - r)^i D_0 \ \text{ for } i \geq 0 -``` - -Equation {eq}`geomseries` expresses $D_i$ as the $i$ th term in the -product of $D_0$ and the geometric series - -$$ -1, (1-r), (1-r)^2, \cdots -$$ - -Therefore, the sum of all deposits in our banking system -$i=0, 1, 2, \ldots$ is - -```{math} -:label: sumdeposits - -\sum_{i=0}^\infty (1-r)^i D_0 = \frac{D_0}{1 - (1-r)} = \frac{D_0}{r} -``` - -### Money Multiplier - -The **money multiplier** is a number that tells the multiplicative -factor by which an exogenous injection of cash into bank $0$ leads -to an increase in the total deposits in the banking system. - -Equation {eq}`sumdeposits` asserts that the **money multiplier** is -$\frac{1}{r}$ - -- An initial deposit of cash of $D_0$ in bank $0$ leads - the banking system to create total deposits of $\frac{D_0}{r}$. -- The initial deposit $D_0$ is held as reserves, distributed - throughout the banking system according to $D_0 = \sum_{i=0}^\infty R_i$. - -## Example: The Keynesian Multiplier - -The famous economist John Maynard Keynes and his followers created a -simple model intended to determine national income $y$ in -circumstances in which - -- there are substantial unemployed resources, in particular **excess - supply** of labor and capital -- prices and interest rates fail to adjust to make aggregate **supply - equal demand** (e.g., prices and interest rates are frozen) -- national income is entirely determined by aggregate demand - -### Static Version - -An elementary Keynesian model of national income determination consists -of three equations that describe aggregate demand for $y$ and its -components. - -The first equation is a national income identity asserting that -consumption $c$ plus investment $i$ equals national income -$y$: - -$$ -c+ i = y -$$ - -The second equation is a Keynesian consumption function asserting that -people consume a fraction $b \in (0,1)$ of their income: - -$$ -c = b y -$$ - -The fraction $b \in (0,1)$ is called the **marginal propensity to -consume**. - -The fraction $1-b \in (0,1)$ is called the **marginal propensity -to save**. - -The third equation simply states that investment is exogenous at level -$i$. - -- *exogenous* means *determined outside this model*. - -Substituting the second equation into the first gives $(1-b) y = i$. - -Solving this equation for $y$ gives - -$$ -y = \frac{1}{1-b} i -$$ - -The quantity $\frac{1}{1-b}$ is called the **investment -multiplier** or simply the **multiplier**. - -Applying the formula for the sum of an infinite geometric series, we can -write the above equation as - -$$ -y = i \sum_{t=0}^\infty b^t -$$ - -where $t$ is a nonnegative integer. - -So we arrive at the following equivalent expressions for the multiplier: - -$$ -\frac{1}{1-b} = \sum_{t=0}^\infty b^t -$$ - -The expression $\sum_{t=0}^\infty b^t$ motivates an interpretation -of the multiplier as the outcome of a dynamic process that we describe -next. - -### Dynamic Version - -We arrive at a dynamic version by interpreting the nonnegative integer -$t$ as indexing time and changing our specification of the -consumption function to take time into account - -- we add a one-period lag in how income affects consumption - -We let $c_t$ be consumption at time $t$ and $i_t$ be -investment at time $t$. - -We modify our consumption function to assume the form - -$$ -c_t = b y_{t-1} -$$ - -so that $b$ is the marginal propensity to consume (now) out of -last period's income. - -We begin with an initial condition stating that - -$$ -y_{-1} = 0 -$$ - -We also assume that - -$$ -i_t = i \ \ \textrm {for all } t \geq 0 -$$ - -so that investment is constant over time. - -It follows that - -$$ -y_0 = i + c_0 = i + b y_{-1} = i -$$ - -and - -$$ -y_1 = c_1 + i = b y_0 + i = (1 + b) i -$$ - -and - -$$ -y_2 = c_2 + i = b y_1 + i = (1 + b + b^2) i -$$ - -and more generally - -$$ -y_t = b y_{t-1} + i = (1+ b + b^2 + \cdots + b^t) i -$$ - -or - -$$ -y_t = \frac{1-b^{t+1}}{1 -b } i -$$ - -Evidently, as $t \rightarrow + \infty$, - -$$ -y_t \rightarrow \frac{1}{1-b} i -$$ - -**Remark 1:** The above formula is often applied to assert that an -exogenous increase in investment of $\Delta i$ at time $0$ -ignites a dynamic process of increases in national income by successive amounts - -$$ -\Delta i, (1 + b )\Delta i, (1+b + b^2) \Delta i , \cdots -$$ - -at times $0, 1, 2, \ldots$. - -**Remark 2** Let $g_t$ be an exogenous sequence of government -expenditures. - -If we generalize the model so that the national income identity -becomes - -$$ -c_t + i_t + g_t = y_t -$$ - -then a version of the preceding argument shows that the **government -expenditures multiplier** is also $\frac{1}{1-b}$, so that a -permanent increase in government expenditures ultimately leads to an -increase in national income equal to the multiplier times the increase -in government expenditures. - -## Example: Interest Rates and Present Values - -We can apply our formula for geometric series to study how interest -rates affect values of streams of dollar payments that extend over time. - -We work in discrete time and assume that $t = 0, 1, 2, \ldots$ -indexes time. - -We let $r \in (0,1)$ be a one-period **net nominal interest rate** - -- if the nominal interest rate is $5$ percent, - then $r= .05$ - -A one-period **gross nominal interest rate** $R$ is defined as - -$$ -R = 1 + r \in (1, 2) -$$ - -- if $r=.05$, then $R = 1.05$ - -**Remark:** The gross nominal interest rate $R$ is an **exchange -rate** or **relative price** of dollars at between times $t$ and -$t+1$. The units of $R$ are dollars at time $t+1$ per -dollar at time $t$. - -When people borrow and lend, they trade dollars now for dollars later or -dollars later for dollars now. - -The price at which these exchanges occur is the gross nominal interest -rate. - -- If I sell $x$ dollars to you today, you pay me $R x$ - dollars tomorrow. -- This means that you borrowed $x$ dollars for me at a gross - interest rate $R$ and a net interest rate $r$. - -We assume that the net nominal interest rate $r$ is fixed over -time, so that $R$ is the gross nominal interest rate at times -$t=0, 1, 2, \ldots$. - -Two important geometric sequences are - -```{math} -:label: geom1 - -1, R, R^2, \cdots -``` - -and - -```{math} -:label: geom2 - -1, R^{-1}, R^{-2}, \cdots -``` - -Sequence {eq}`geom1` tells us how dollar values of an investment **accumulate** -through time. - -Sequence {eq}`geom2` tells us how to **discount** future dollars to get their -values in terms of today's dollars. - -### Accumulation - -Geometric sequence {eq}`geom1` tells us how one dollar invested and re-invested -in a project with gross one period nominal rate of return accumulates - -- here we assume that net interest payments are reinvested in the - project -- thus, $1$ dollar invested at time $0$ pays interest - $r$ dollars after one period, so we have $r+1 = R$ - dollars at time$1$ -- at time $1$ we reinvest $1+r =R$ dollars and receive interest - of $r R$ dollars at time $2$ plus the *principal* - $R$ dollars, so we receive $r R + R = (1+r)R = R^2$ - dollars at the end of period $2$ -- and so on - -Evidently, if we invest $x$ dollars at time $0$ and -reinvest the proceeds, then the sequence - -$$ -x , xR , x R^2, \cdots -$$ - -tells how our account accumulates at dates $t=0, 1, 2, \ldots$. - -### Discounting - -Geometric sequence {eq}`geom2` tells us how much future dollars are worth in terms of today's dollars. - -Remember that the units of $R$ are dollars at $t+1$ per -dollar at $t$. - -It follows that - -- the units of $R^{-1}$ are dollars at $t$ per dollar at $t+1$ -- the units of $R^{-2}$ are dollars at $t$ per dollar at $t+2$ -- and so on; the units of $R^{-j}$ are dollars at $t$ per - dollar at $t+j$ - -So if someone has a claim on $x$ dollars at time $t+j$, it -is worth $x R^{-j}$ dollars at time $t$ (e.g., today). - -### Application to Asset Pricing - -A **lease** requires a payments stream of $x_t$ dollars at -times $t = 0, 1, 2, \ldots$ where - -$$ -x_t = G^t x_0 -$$ - -where $G = (1+g)$ and $g \in (0,1)$. - -Thus, lease payments increase at $g$ percent per period. - -For a reason soon to be revealed, we assume that $G < R$. - -The **present value** of the lease is - -$$ -\begin{aligned} p_0 & = x_0 + x_1/R + x_2/(R^2) + \ddots \\ - & = x_0 (1 + G R^{-1} + G^2 R^{-2} + \cdots ) \\ - & = x_0 \frac{1}{1 - G R^{-1}} \end{aligned} -$$ - -where the last line uses the formula for an infinite geometric series. - -Recall that $R = 1+r$ and $G = 1+g$ and that $R > G$ -and $r > g$ and that $r$ and $g$ are typically small -numbers, e.g., .05 or .03. - -Use the Taylor series of $\frac{1}{1+r}$ about $r=0$, -namely, - -$$ -\frac{1}{1+r} = 1 - r + r^2 - r^3 + \cdots -$$ - -and the fact that $r$ is small to approximate -$\frac{1}{1+r} \approx 1 - r$. - -Use this approximation to write $p_0$ as - -$$ -\begin{aligned} - p_0 &= x_0 \frac{1}{1 - G R^{-1}} \\ - &= x_0 \frac{1}{1 - (1+g) (1-r) } \\ - &= x_0 \frac{1}{1 - (1+g - r - rg)} \\ - & \approx x_0 \frac{1}{r -g } -\end{aligned} -$$ - -where the last step uses the approximation $r g \approx 0$. - -The approximation - -$$ -p_0 = \frac{x_0 }{r -g } -$$ - -is known as the **Gordon formula** for the present value or current -price of an infinite payment stream $x_0 G^t$ when the nominal -one-period interest rate is $r$ and when $r > g$. - -We can also extend the asset pricing formula so that it applies to finite leases. - -Let the payment stream on the lease now be $x_t$ for $t= 1,2, \dots,T$, where again - -$$ -x_t = G^t x_0 -$$ - -The present value of this lease is: - -$$ -\begin{aligned} \begin{split}p_0&=x_0 + x_1/R + \dots +x_T/R^T \\ &= x_0(1+GR^{-1}+\dots +G^{T}R^{-T}) \\ &= \frac{x_0(1-G^{T+1}R^{-(T+1)})}{1-GR^{-1}} \end{split}\end{aligned} -$$ - -Applying the Taylor series to $R^{-(T+1)}$ about $r=0$ we get: - -$$ -\frac{1}{(1+r)^{T+1}}= 1-r(T+1)+\frac{1}{2}r^2(T+1)(T+2)+\dots \approx 1-r(T+1) -$$ - -Similarly, applying the Taylor series to $G^{T+1}$ about $g=0$: - -$$ -(1+g)^{T+1} = 1+(T+1)g+\frac{T(T+1)}{2!}g^2+\frac{(T-1)T(T+1)}{3!}g^3+\dots \approx 1+ (T+1)g -$$ - -Thus, we get the following approximation: - -$$ -p_0 =\frac{x_0(1-(1+(T+1)g)(1-r(T+1)))}{1-(1-r)(1+g) } -$$ - -Expanding: - -$$ -\begin{aligned} p_0 &=\frac{x_0(1-1+(T+1)^2 rg -r(T+1)+g(T+1))}{1-1+r-g+rg} \\&=\frac{x_0(T+1)((T+1)rg+r-g)}{r-g+rg} \\ &\approx \frac{x_0(T+1)(r-g)}{r-g}+\frac{x_0rg(T+1)}{r-g}\\ &= x_0(T+1) + \frac{x_0rg(T+1)}{r-g} \end{aligned} -$$ - -We could have also approximated by removing the second term -$rgx_0(T+1)$ when $T$ is relatively small compared to -$1/(rg)$ to get $x_0(T+1)$ as in the finite stream -approximation. - -We will plot the true finite stream present-value and the two -approximations, under different values of $T$, and $g$ and $r$ in Python. - -First we plot the true finite stream present-value after computing it -below - -```{code-cell} python3 -# True present value of a finite lease -def finite_lease_pv_true(T, g, r, x_0): - G = (1 + g) - R = (1 + r) - return (x_0 * (1 - G**(T + 1) * R**(-T - 1))) / (1 - G * R**(-1)) -# First approximation for our finite lease - -def finite_lease_pv_approx_1(T, g, r, x_0): - p = x_0 * (T + 1) + x_0 * r * g * (T + 1) / (r - g) - return p - -# Second approximation for our finite lease -def finite_lease_pv_approx_2(T, g, r, x_0): - return (x_0 * (T + 1)) - -# Infinite lease -def infinite_lease(g, r, x_0): - G = (1 + g) - R = (1 + r) - return x_0 / (1 - G * R**(-1)) -``` - -Now that we have defined our functions, we can plot some outcomes. - -First we study the quality of our approximations - -```{code-cell} python3 -def plot_function(axes, x_vals, func, args): - axes.plot(x_vals, func(*args), label=func.__name__) - -T_max = 50 - -T = np.arange(0, T_max+1) -g = 0.02 -r = 0.03 -x_0 = 1 - -our_args = (T, g, r, x_0) -funcs = [finite_lease_pv_true, - finite_lease_pv_approx_1, - finite_lease_pv_approx_2] - ## the three functions we want to compare - -fig, ax = plt.subplots() -ax.set_title('Finite Lease Present Value $T$ Periods Ahead') -for f in funcs: - plot_function(ax, T, f, our_args) -ax.legend() -ax.set_xlabel('$T$ Periods Ahead') -ax.set_ylabel('Present Value, $p_0$') -plt.show() -``` - -Evidently our approximations perform well for small values of $T$. - -However, holding $g$ and r fixed, our approximations deteriorate as $T$ increases. - -Next we compare the infinite and finite duration lease present values -over different lease lengths $T$. - -```{code-cell} python3 -# Convergence of infinite and finite -T_max = 1000 -T = np.arange(0, T_max+1) -fig, ax = plt.subplots() -ax.set_title('Infinite and Finite Lease Present Value $T$ Periods Ahead') -f_1 = finite_lease_pv_true(T, g, r, x_0) -f_2 = np.full(T_max+1, infinite_lease(g, r, x_0)) -ax.plot(T, f_1, label='T-period lease PV') -ax.plot(T, f_2, '--', label='Infinite lease PV') -ax.set_xlabel('$T$ Periods Ahead') -ax.set_ylabel('Present Value, $p_0$') -ax.legend() -plt.show() -``` - -The graph above shows how as duration $T \rightarrow +\infty$, -the value of a lease of duration $T$ approaches the value of a -perpetual lease. - -Now we consider two different views of what happens as $r$ and -$g$ covary - -```{code-cell} python3 -# First view -# Changing r and g -fig, ax = plt.subplots() -ax.set_title('Value of lease of length $T$') -ax.set_ylabel('Present Value, $p_0$') -ax.set_xlabel('$T$ periods ahead') -T_max = 10 -T=np.arange(0, T_max+1) - -rs, gs = (0.9, 0.5, 0.4001, 0.4), (0.4, 0.4, 0.4, 0.5), -comparisons = ('$\gg$', '$>$', r'$\approx$', '$<$') -for r, g, comp in zip(rs, gs, comparisons): - ax.plot(finite_lease_pv_true(T, g, r, x_0), label=f'r(={r}) {comp} g(={g})') - -ax.legend() -plt.show() -``` - -This graph gives a big hint for why the condition $r > g$ is -necessary if a lease of length $T = +\infty$ is to have finite -value. - -For fans of 3-d graphs the same point comes through in the following -graph. - -If you aren't enamored of 3-d graphs, feel free to skip the next -visualization! - -```{code-cell} python3 -# Second view -fig = plt.figure() -T = 3 -ax = plt.subplot(projection='3d') -r = np.arange(0.01, 0.99, 0.005) -g = np.arange(0.011, 0.991, 0.005) - -rr, gg = np.meshgrid(r, g) -z = finite_lease_pv_true(T, gg, rr, x_0) - -# Removes points where undefined -same = (rr == gg) -z[same] = np.nan -surf = ax.plot_surface(rr, gg, z, cmap=cm.coolwarm, - antialiased=True, clim=(0, 15)) -fig.colorbar(surf, shrink=0.5, aspect=5) -ax.set_xlabel('$r$') -ax.set_ylabel('$g$') -ax.set_zlabel('Present Value, $p_0$') -ax.view_init(20, 10) -ax.set_title('Three Period Lease PV with Varying $g$ and $r$') -plt.show() -``` - -We can use a little calculus to study how the present value $p_0$ -of a lease varies with $r$ and $g$. - -We will use a library called [SymPy](https://www.sympy.org/). - -SymPy enables us to do symbolic math calculations including -computing derivatives of algebraic equations. - -We will illustrate how it works by creating a symbolic expression that -represents our present value formula for an infinite lease. - -After that, we'll use SymPy to compute derivatives - -```{code-cell} python3 -# Creates algebraic symbols that can be used in an algebraic expression -g, r, x0 = sym.symbols('g, r, x0') -G = (1 + g) -R = (1 + r) -p0 = x0 / (1 - G * R**(-1)) -init_printing(use_latex='mathjax') -print('Our formula is:') -p0 -``` - -```{code-cell} python3 -print('dp0 / dg is:') -dp_dg = sym.diff(p0, g) -dp_dg -``` - -```{code-cell} python3 -print('dp0 / dr is:') -dp_dr = sym.diff(p0, r) -dp_dr -``` - -We can see that for $\frac{\partial p_0}{\partial r}<0$ as long as -$r>g$, $r>0$ and $g>0$ and $x_0$ is positive, -so $\frac{\partial p_0}{\partial r}$ will always be negative. - -Similarly, $\frac{\partial p_0}{\partial g}>0$ as long as $r>g$, $r>0$ and $g>0$ and $x_0$ is positive, so $\frac{\partial p_0}{\partial g}$ -will always be positive. - -## Back to the Keynesian Multiplier - -We will now go back to the case of the Keynesian multiplier and plot the -time path of $y_t$, given that consumption is a constant fraction -of national income, and investment is fixed. - -```{code-cell} python3 -# Function that calculates a path of y -def calculate_y(i, b, g, T, y_init): - y = np.zeros(T+1) - y[0] = i + b * y_init + g - for t in range(1, T+1): - y[t] = b * y[t-1] + i + g - return y - -# Initial values -i_0 = 0.3 -g_0 = 0.3 -# 2/3 of income goes towards consumption -b = 2/3 -y_init = 0 -T = 100 - -fig, ax = plt.subplots() -ax.set_title('Path of Aggregate Output Over Time') -ax.set_xlabel('$t$') -ax.set_ylabel('$y_t$') -ax.plot(np.arange(0, T+1), calculate_y(i_0, b, g_0, T, y_init)) -# Output predicted by geometric series -ax.hlines(i_0 / (1 - b) + g_0 / (1 - b), xmin=-1, xmax=101, linestyles='--') -plt.show() -``` - -In this model, income grows over time, until it gradually converges to -the infinite geometric series sum of income. - -We now examine what will -happen if we vary the so-called **marginal propensity to consume**, -i.e., the fraction of income that is consumed - -```{code-cell} python3 -bs = (1/3, 2/3, 5/6, 0.9) - -fig,ax = plt.subplots() -ax.set_title('Changing Consumption as a Fraction of Income') -ax.set_ylabel('$y_t$') -ax.set_xlabel('$t$') -x = np.arange(0, T+1) -for b in bs: - y = calculate_y(i_0, b, g_0, T, y_init) - ax.plot(x, y, label=r'$b=$'+f"{b:.2f}") -ax.legend() -plt.show() -``` - -Increasing the marginal propensity to consume $b$ increases the -path of output over time. - -Now we will compare the effects on output of increases in investment and government spending. - -```{code-cell} python3 -fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 10)) -fig.subplots_adjust(hspace=0.3) - -x = np.arange(0, T+1) -values = [0.3, 0.4] - -for i in values: - y = calculate_y(i, b, g_0, T, y_init) - ax1.plot(x, y, label=f"i={i}") -for g in values: - y = calculate_y(i_0, b, g, T, y_init) - ax2.plot(x, y, label=f"g={g}") - -axes = ax1, ax2 -param_labels = "Investment", "Government Spending" -for ax, param in zip(axes, param_labels): - ax.set_title(f'An Increase in {param} on Output') - ax.legend(loc ="lower right") - ax.set_ylabel('$y_t$') - ax.set_xlabel('$t$') -plt.show() -``` - -Notice here, whether government spending increases from 0.3 to 0.4 or -investment increases from 0.3 to 0.4, the shifts in the graphs are -identical. diff --git a/lectures/lp_intro.md b/lectures/lp_intro.md deleted file mode 100644 index 39e551218..000000000 --- a/lectures/lp_intro.md +++ /dev/null @@ -1,808 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.14.4 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# Linear Programming - -## Overview - -**Linear programming** problems either maximize or minimize -a linear objective function subject to a set of linear equality and/or inequality constraints. - -Linear programs come in pairs: - -* an original **primal** problem, and - -* an associated **dual** problem. - -If a primal problem involves **maximization**, the dual problem involves **minimization**. - -If a primal problem involves **minimization**, the dual problem involves **maximization**. - -We provide a standard form of a linear program and methods to transform other forms of linear programming problems into a standard form. - -We tell how to solve a linear programming problem using [SciPy](https://scipy.org/). - -We describe the important concept of complementary slackness and how it relates to the dual problem. - -Let's start with some standard imports. - -```{code-cell} ipython3 -import numpy as np -from scipy.optimize import linprog -import matplotlib.pyplot as plt -from matplotlib.patches import Polygon -``` - -## Objective Function and Constraints - -We want to minimize a **cost function** $c'x = \sum_{i=1}^n c_i x_i$ over feasible values of $x = (x_1,x_2,\dots,x_n)'$. - -Here - -* $c = (c_1,c_2,\dots,c_n)'$ is a **unit cost vector**, and - -* $x = (x_1,x_2,\dots,x_n)'$ is a vector of **decision variables** - -Decision variables are restricted to satisfy a set of linear equality and/or inequality constraints. - -We describe the constraints with the following collections of $n$-dimensional vectors $a_i$ and scalars $b_i$ and associated sets indexing the equality and inequality constraints: - -* $a_i$ for $i \in M_i$, where $M_1,M_2,M_3$ are each sets of indexes - -and a collection of scalers - -* $b_i$ for $i \in N_i$, where $N_1,N_2,N_3$ are each sets of indexes. - -A linear programming can be stated as {cite}`bertsimas_tsitsiklis1997`: - -$$ -\begin{aligned} -\min_{x} \ & c' x \\ -\mbox{subject to } \ & a_i' x \ge b_i, & i \in M_1 \\ -& a_i' x \le b_i, & i \in M_2 \\ -& a_i' x = b_i, & i \in M_3 \\ -& x_j \ge 0, & j \in N_1 \\ -& x_j \le 0, & j \in N_2 \\ -& x_j\ \text{unrestricted}, & j \in N_3 \\ -\end{aligned} -$$ (linprog) - -A vector $x$ that satisfies all of the constraints is called a **feasible solution**. - -A collection of all feasible solutions is called a **feasible set**. - -A feasible solution $x$ that minimizes the cost function is called an **optimal solution**. - -The corresponding value of cost function $c'x$ is called the **optimal value**. - -If the feasible set is empty, we say that solving the linear programming problem is **infeasible**. - -If, for any $K \in \mathbb R$, there exists a feasible solution $x$ such that $c'x < K$, we say that the problem is **unbounded** or equivalently that the optimal value is $-\infty$. - -## Example 1: Production Problem - -This example was created by {cite}`bertsimas_tsitsiklis1997` - -Suppose that a factory can produce two goods called Product $1$ and Product $2$. - -To produce each product requires both material and labor. - -Selling each product generates revenue. - -Required per unit material and labor inputs and revenues are shown in table below: - -| | Product 1 | Product 2 | -| :------: | :-------: | :-------: | -| Material | 2 | 5 | -| Labor | 4 | 2 | -| Revenue | 3 | 4 | - -30 units of material and 20 units of labor available. - -A firm's problem is to construct a production plan that uses its 30 units of materials and 20 unites of labor -to maximize its revenue. - -Let $x_i$ denote the quantity of Product $i$ that the firm produces. - -This problem can be formulated as: - -$$ -\begin{aligned} -\max_{x_1,x_2} \ & z = 3 x_1 + 4 x_2 \\ -\mbox{subject to } \ & 2 x_1 + 5 x_2 \le 30 \\ -& 4 x_1 + 2 x_2 \le 20 \\ -& x_1, x_2 \ge 0 \\ -\end{aligned} -$$ - -The following graph illustrates the firm's constraints and iso-revenue lines. - -```{code-cell} ipython3 -fig, ax = plt.subplots(figsize=(8, 6)) -ax.grid() - -# Draw constraint lines -ax.hlines(0, -1, 17.5) -ax.vlines(0, -1, 12) -ax.plot(np.linspace(-1, 17.5, 100), 6-0.4*np.linspace(-1, 17.5, 100), color="c") -ax.plot(np.linspace(-1, 5.5, 100), 10-2*np.linspace(-1, 5.5, 100), color="c") -ax.text(1.5, 8, "$2x_1 + 5x_2 \leq 30$", size=12) -ax.text(10, 2.5, "$4x_1 + 2x_2 \leq 20$", size=12) -ax.text(-2, 2, "$x_2 \geq 0$", size=12) -ax.text(2.5, -0.7, "$x_1 \geq 0$", size=12) - -# Draw the feasible region -feasible_set = Polygon(np.array([[0, 0], - [0, 6], - [2.5, 5], - [5, 0]]), - color="cyan") -ax.add_patch(feasible_set) - -# Draw the objective function -ax.plot(np.linspace(-1, 5.5, 100), 3.875-0.75*np.linspace(-1, 5.5, 100), color="orange") -ax.plot(np.linspace(-1, 5.5, 100), 5.375-0.75*np.linspace(-1, 5.5, 100), color="orange") -ax.plot(np.linspace(-1, 5.5, 100), 6.875-0.75*np.linspace(-1, 5.5, 100), color="orange") -ax.arrow(-1.6, 5, 0, 2, width = 0.05, head_width=0.2, head_length=0.5, color="orange") -ax.text(5.7, 1, "$z = 3x_1 + 4x_2$", size=12) - -# Draw the optimal solution -ax.plot(2.5, 5, "*", color="black") -ax.text(2.7, 5.2, "Optimal Solution", size=12) - -plt.show() -``` - -The blue region is the feasible set within which all constraints are satisfied. - -Parallel orange lines are iso-revenue lines. - -The firm's objective is to find the parallel orange lines to the upper boundary of the feasible set. - -The intersection of the feasible set and the highest orange line delineates the optimal set. - -In this example, the optimal set is the point $(2.5, 5)$. - -## Example 2: Investment Problem - -We now consider a problem posed and solved by {cite}`hu_guo2018`. - -A mutual fund has $ \$ 100,000$ to be invested over a three year horizon. - -Three investment options are available: - -1. **Annuity:** the fund can pay a same amount of new capital at the beginning of each of three years and receive a payoff of 130\% of **total capital** invested at the end of the third year. Once the mutual fund decides to invest in this annuity, it has to keep investing in all subsequent years in the three year horizon. - -2. **Bank account:** the fund can deposit any amount into a bank at the beginning of each year and receive its capital plus 6\% interest at the end of that year. In addition, the mutual fund is permitted to borrow no more than $20,000 at the beginning of each year and is asked to pay back the amount borrowed plus 6\% interest at the end of the year. The mutual fund can choose whether to deposit or borrow at the beginning of each year. - -3. **Corporate bond:** At the beginning of the second year, a corporate bond becomes available. -The fund can buy an amount -that is no more than $ \$ $50,000 of this bond at the beginning of the second year and at the end of the third year receive a payout of 130\% of the amount invested in the bond. - -The mutual fund's objective is to maximize total payout that it owns at the end of the third year. - -We can formulate this as a linear programming problem. - -Let $x_1$ be the amount of put in the annuity, $x_2, x_3, x_4$ be bank deposit balances at the beginning of the three years, and $x_5$ be the amount invested in the corporate bond. - -When $x_2, x_3, x_4$ are negative, it means that the mutual fund has borrowed from bank. - -The table below shows the mutual fund's decision variables together with the timing protocol described above: - -| | Year 1 | Year 2 | Year 3 | -| :------------: | :----: | :----: | :----: | -| Annuity | $x_1$ | $x_1$ | $x_1$ | -| Bank account | $x_2$ | $x_3$ | $x_4$ | -| Corporate bond | 0 | $x_5$ | 0 | - -The mutual fund's decision making proceeds according to the following timing protocol: - -1. At the beginning of the first year, the mutual fund decides how much to invest in the annuity and - how much to deposit in the bank. This decision is subject to the constraint: - - $$ - x_1 + x_2 = 100,000 - $$ - -2. At the beginning of the second year, the mutual fund has a bank balance of $1.06 x_2$. - It must keep $x_1$ in the annuity. It can choose to put $x_5$ into the corporate bond, - and put $x_3$ in the bank. These decisions are restricted by - - $$ - x_1 + x_5 = 1.06 x_2 - x_3 - $$ - -3. At the beginning of the third year, the mutual fund has a bank account balance equal - to $1.06 x_3$. It must again invest $x_1$ in the annuity, - leaving it with a bank account balance equal to $x_4$. This situation is summarized by the restriction: - - $$ - x_1 = 1.06 x_3 - x_4 - $$ - -The mutual fund's objective function, i.e., its wealth at the end of the third year is: - -$$ -1.30 \cdot 3x_1 + 1.06 x_4 + 1.30 x_5 -$$ - -Thus, the mutual fund confronts the linear program: - -$$ -\begin{aligned} -\max_{x} \ & 1.30 \cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\ -\mbox{subject to } \ & x_1 + x_2 = 100,000\\ - & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\ - & x_1 - 1.06 x_3 + x_4 = 0\\ - & x_2 \ge -20,000\\ - & x_3 \ge -20,000\\ - & x_4 \ge -20,000\\ - & x_5 \le 50,000\\ - & x_j \ge 0, \quad j = 1,5\\ - & x_j \ \text{unrestricted}, \quad j = 2,3,4\\ -\end{aligned} -$$ - -## Standard Form - -For purposes of - -* unifying linear programs that are initially stated in superficially different forms, and - -* having a form that is convenient to put into black-box software packages, - -it is useful to devote some effort to describe a **standard form**. - -Our standard form is: - -$$ -\begin{aligned} -\min_{x} \ & c_1 x_1 + c_2 x_2 + \dots + c_n x_n \\ -\mbox{subject to } \ & a_{11} x_1 + a_{12} x_2 + \dots + a_{1n} x_n = b_1 \\ - & a_{21} x_1 + a_{22} x_2 + \dots + a_{2n} x_n = b_2 \\ - & \quad \vdots \\ - & a_{m1} x_1 + a_{m2} x_2 + \dots + a_{mn} x_n = b_m \\ - & x_1, x_2, \dots, x_n \ge 0 \\ -\end{aligned} -$$ - -Let - -$$ -A = \begin{bmatrix} -a_{11} & a_{12} & \dots & a_{1n} \\ -a_{21} & a_{22} & \dots & a_{2n} \\ - & & \vdots & \\ -a_{m1} & a_{m2} & \dots & a_{mn} \\ -\end{bmatrix}, \quad -b = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \\ \end{bmatrix}, \quad -c = \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \\ \end{bmatrix}, \quad -x = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \\ \end{bmatrix}. \quad -$$ - -The standard form LP problem can be expressed concisely as: - -$$ -\begin{aligned} -\min_{x} \ & c'x \\ -\mbox{subject to } \ & Ax = b\\ - & x >= 0\\ -\end{aligned} -$$ (lpproblem) - -Here, $Ax = b$ means that the $i$-th entry of $Ax$ equals the $i$-th entry of $b$ for every $i$. - -Similarly, $x >= 0$ means that $x_j$ is greater than $0$ for every $j$. - -### Useful Transformations - -It is useful to know how to transform a problem that initially is not stated in the standard form into one that is. - -By deploying the following steps, any linear programming problem can be transformed into an equivalent standard form linear programming problem. - -1. **Objective Function:** If a problem is originally a constrained **maximization** problem, we can construct a new objective function that is the additive inverse of the original objective function. The transformed problem is then a **minimization** problem. - -2. **Decision Variables:** Given a variable $x_j$ satisfying $x_j \le 0$, we can introduce a new variable $x_j' = - x_j$ and subsitute it into original problem. Given a free variable $x_i$ with no restriction on its sign, we can introduce two new variables $x_j^+$ and $x_j^-$ satisfying $x_j^+, x_j^- \ge 0$ and replace $x_j$ by $x_j^+ - x_j^-$. - -3. **Inequality constraints:** Given an inequality constraint $\sum_{j=1}^n a_{ij}x_j \le 0$, we can introduce a new variable $s_i$, called a **slack variable** that satisfies $s_i \ge 0$ and replace the original constraint by $\sum_{j=1}^n a_{ij}x_j + s_i = 0$. - -Let's apply the above steps to the two examples described above. - -### Example 1: Production Problem - -The original problem is: - -$$ -\begin{aligned} -\max_{x_1,x_2} \ & 3 x_1 + 4 x_2 \\ -\mbox{subject to } \ & 2 x_1 + 5 x_2 \le 30 \\ -& 4 x_1 + 2 x_2 \le 20 \\ -& x_1, x_2 \ge 0 \\ -\end{aligned} -$$ - -This problem is equivalent to the following problem with a standard form: - -$$ -\begin{aligned} -\min_{x_1,x_2} \ & -(3 x_1 + 4 x_2) \\ -\mbox{subject to } \ & 2 x_1 + 5 x_2 + s_1 = 30 \\ -& 4 x_1 + 2 x_2 + s_2 = 20 \\ -& x_1, x_2, s_1, s_2 \ge 0 \\ -\end{aligned} -$$ - -### Example 2: Investment Problem - -The original problem is: - -$$ -\begin{aligned} -\max_{x} \ & 1.30 \cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\ -\mbox{subject to } \ & x_1 + x_2 = 100,000\\ - & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\ - & x_1 - 1.06 x_3 + x_4 = 0\\ - & x_2 \ge -20,000\\ - & x_3 \ge -20,000\\ - & x_4 \ge -20,000\\ - & x_5 \le 50,000\\ - & x_j \ge 0, \quad j = 1,5\\ - & x_j \ \text{unrestricted}, \quad j = 2,3,4\\ -\end{aligned} -$$ - -This problem is equivalent to the following problem with a standard form: - -$$ -\begin{aligned} -\min_{x} \ & -(1.30 \cdot 3x_1 + 1.06 x_4^+ - 1.06 x_4^- + 1.30 x_5) \\ -\mbox{subject to } \ & x_1 + x_2^+ - x_2^- = 100,000\\ - & x_1 - 1.06 (x_2^+ - x_2^-) + x_3^+ - x_3^- + x_5 = 0\\ - & x_1 - 1.06 (x_3^+ - x_3^-) + x_4^+ - x_4^- = 0\\ - & x_2^- - x_2^+ + s_1 = 20,000\\ - & x_3^- - x_3^+ + s_2 = 20,000\\ - & x_4^- - x_4^+ + s_3 = 20,000\\ - & x_5 + s_4 = 50,000\\ - & x_j \ge 0, \quad j = 1,5\\ - & x_j^+, x_j^- \ge 0, \quad j = 2,3,4\\ - & s_j \ge 0, \quad j = 1,2,3,4\\ -\end{aligned} -$$ - -## Computations - -The package *scipy.optimize* provides a function ***linprog*** to solve linear programming problems with a form below: - -$$ -\begin{aligned} -\min_{x} \ & c' x \\ -\mbox{subject to } \ & A_{ub}x \le b_{ub} \\ - & A_{eq}x = b_{eq} \\ - & l \le x \le u \\ -\end{aligned} -$$ - -```{note} -By default $l = 0$ and $u = \text{None}$ unless explicitly specified with the argument 'bounds'. -``` - -Let's apply this great Python tool to solve our two example problems. - -### Example 1: Production Problem - -The problem is: - -$$ -\begin{aligned} -\max_{x_1,x_2} \ & 3 x_1 + 4 x_2 \\ -\mbox{subject to } \ & 2 x_1 + 5 x_2 \le 30 \\ -& 4 x_1 + 2 x_2 \le 20 \\ -& x_1, x_2 \ge 0 \\ -\end{aligned} -$$ - -```{code-cell} ipython3 -# Construct parameters -c_ex1 = np.array([3, 4]) - -# Inequality constraints -A_ex1 = np.array([[2, 5], - [4, 2]]) -b_ex1 = np.array([30,20]) - -# Solve the problem -# we put a negative sign on the objective as linprog does minimization -res_ex1 = linprog(-c_ex1, A_ub=A_ex1, b_ub=b_ex1) - -res_ex1 -``` - -The optimal plan tells the factory to produce 2.5 units of Product 1 and 5 units of Product 2; that generates a maximizing value of revenue of 27.5. - -We are using the *linprog* function as a **black box**. - -Inside it, Python first transforms the problem into standard form. - -To do that, for each inequality constraint it generates one slack variable. - -Here the vector of slack variables is a two-dimensional NumPy array that equals $b_{ub} - A_{ub}x$. - -See the [official documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html#scipy.optimize.linprog) for more details. - -```{note} -This problem is to maximize the objective, so that we need to put a minus sign in front of parameter vector c. -``` - -### Example 2: Investment Problem - -The problem is: - -$$ -\begin{aligned} -\max_{x} \ & 1.30 \cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\ -\mbox{subject to } \ & x_1 + x_2 = 100,000\\ - & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\ - & x_1 - 1.06 x_3 + x_4 = 0\\ - & x_2 \ge -20,000\\ - & x_3 \ge -20,000\\ - & x_4 \ge -20,000\\ - & x_5 \le 50,000\\ - & x_j \ge 0, \quad j = 1,5\\ - & x_j \ \text{unrestricted}, \quad j = 2,3,4\\ -\end{aligned} -$$ - -Let's solve this problem using *linprog*. - -```{code-cell} ipython3 -# Construct parameters -rate = 1.06 - -# Objective function parameters -c_ex2 = np.array([1.30*3, 0, 0, 1.06, 1.30]) - -# Inequality constraints -A_ex2 = np.array([[1, 1, 0, 0, 0], - [1, -rate, 1, 0, 1], - [1, 0, -rate, 1, 0]]) -b_ex2 = np.array([100000, 0, 0]) - -# Bounds on decision variables -bounds_ex2 = [( 0, None), - (-20000, None), - (-20000, None), - (-20000, None), - ( 0, 50000)] - -# Solve the problem -res_ex2 = linprog(-c_ex2, A_eq=A_ex2, b_eq=b_ex2, - bounds=bounds_ex2) - -res_ex2 -``` - -Python tells us that the best investment strategy is: - -1. At the beginning of the first year, the mutual fund should buy $ \$24,927.75$ of the annuity. Its bank account balance should be $ \$75,072.25$. - -2. At the beginning of the second year, the mutual fund should buy $ \$50,000 $ of the corporate bond and keep invest in the annuity. Its bank account balance should be $ \$ 4,648.83$. - -3. At the beginning of the third year, the mutual fund should borrow $ \$20,000$ from the bank and invest in the annuity. - -4. At the end of the third year, the mutual fund will get payouts from the annuity and corporate bond and repay its loan from the bank. At the end it will own $ \$ $141018.24, so that it's total net rate of return over the three periods is 41.02\%. - -## Duality - -Associated with a linear programming of form {eq}`linprog` with $m$ constraints and $n$ decision variables, -there is an **dual** linear programming problem that takes the form (please see {cite}`bertsimas_tsitsiklis1997`) - -$$ -\begin{aligned} -\max_{p} \ & b' p \\ -\mbox{subject to } \ & p_i \ge 0, & i \in M_1 \\ -& p_i \le 0, & i \in M_2 \\ -& p_i\ \text{unrestricted}, & i \in M_3 \\ -& A_j' p \le c_j, & j \in N_1 \\ -& A_j' p \ge c_j, & j \in N_2 \\ -& A_j' p = c_j, & j \in N_3 \\ -\end{aligned} -$$ - -Where $A_j$ is $j$-th column of the $m$ by $n$ matrix $A$. - -```{note} -In what follows, we shall use $a_i'$ to denote the $i$-th row of $A$ and $A_j$ to denote the $j$-th column of $A$. -``` - -$$ -A = \begin{bmatrix} -a_1' \\ -a_2' \\ -\ \\ -a_m' \\ -\end{bmatrix}. -$$ - -To construct the dual of linear programming problem {eq}`linprog`, we proceed as follows: - -1. For every constraint $a_i' x \ge(\le or =) b_i$, $j = 1,2,...,m$, in the primal problem, we construct a corresponding dual variable $p_i$. $p_i$ is restricted to be positive if $a_i' x \ge b_i$ or negative if $a_i' x \le b_i$ or unrestricted if $a_i' x = b_i$. We construct the $m$-dimensional vector $p$ with entries $p_i$. - -2. For every variable $x_j$, $j = 1,2,...,n$, we construct a corresponding dual constraint $A_j' p \ge(\le or =) c_j$. The constraint is $A_j' p \ge c_j$ if $x_j \le 0$, $A_j' p \le c_j$ if $x_j \ge 0$ or $A_j' p = c_j$ if $x_j $ is unrestricted. - -3. The dual problem is to **maximize** objective function $b'p$. - -For a **maximization** problem, we can first transform it to an equivalent minimization problem and then follow the above steps above to construct the dual **minimization** problem. - -We can easily verify that **the dual of a dual problem is the primal problem**. - -The following table summarizes relationships between objects in primal and dual problems. - -| Objective: Min | Objective: Max | -| :--------------: | :--------------: | -| m constraints | m variables | -| constraint $\ge$ | variable $\ge$ 0 | -| constraint $\le$ | variable $\le$ 0 | -| constraint $=$ | variable free | -| n variables | n constraints | -| variable $\ge$ 0 | constraint $\le$ | -| variable $\le$ 0 | constraint $\ge$ | -| variable free | constraint $=$ | - -As an example, the dual problem of the standard form {eq}`lpproblem` is: - -$$ -\begin{aligned} -\max_{p} \ & b'p \\ -\mbox{subject to } \ & A'p \le c\\ -\end{aligned} -$$ - -As another example, consider a linear programming problem with form: - -$$ -\begin{aligned} -\max_{x} \ & c'x \\ -\mbox{subject to } \ & A x \le b\\ -& x \ge 0\\ -\end{aligned} -$$ (linprog2) - -Its dual problem is: - -$$ -\begin{aligned} -\min_{p} \ & b'p \\ -\mbox{subject to } \ & A' p \ge c\\ -& p \ge 0\\ -\end{aligned} -$$ - -## Duality Theorems - -Primal and dual problems are linked by powerful **duality theorems** that have **weak** and **strong** forms. - -The duality theorems provide the foundations of enlightening economic interpretations of linear programming problems. - -**Weak duality:** For linear programming problem {eq}`linprog`, if $x$ and $p$ are feasible solutions to the primal and the dual problems, respectively, then - -$$ -b'p \le c'x -$$ - -**Strong duality:** For linear programming problem {eq}`linprog`, if the primal problem has an optimal solution $x$, then the dual problem also has an optimal solution. Denote an optimal solution of the dual problem as $p$. Then - -$$ -b'p = c'x -$$ - -According to strong duality, we can find the optimal value for the primal problem by solving the dual problem. - -But the dual problem tells us even more as we shall see next. - -### Complementary Slackness - -Let $x$ and $p$ be feasible solutions to the primal problem {eq}`linprog` and its dual problem, respectively. - -Then $x$ and $p$ are also optimal solutions of the primal and dual problems if and only if: - -$$ -p_i (a_i' x - b_i) = 0, \quad \forall i, \\ -x_j (A_j' p - c_j) = 0, \quad \forall j. -$$ - -This means that $p_i = 0$ if $a_i' x - b_i \neq 0$ and $x_j = 0$ if $A_j' p - c_j \neq 0$. - -These are the celebrated **complementary slackness** conditions. - -Let's interpret them. - -### Interpretations - -Let's take a version of problem {eq}`linprog2` as a production problem and consider its associated dual problem. - -A factory produce $n$ products with $m$ types of resources. - -Where $i=1,2,\dots,m$ and $j=1,2,\dots,n$, let - -* $x_j$ denote quantities of product $j$ to be produced - -* $a_{ij}$ denote required amount of resource $i$ to make one unit of product $j$, - -* $b_i$ denotes the avaliable amount of resource $i$ - -* $c_j$ denotes the revenue generated by producing one unit of product $j$. - - -**Dual variables:** By strong duality, we have - -$$ -c_1 x_1 + c_2 x_2 + \dots + c_n x_n = b_1 p_1 + b_2 p_2 + \dots + b_m p_m. -$$ - -Evidently, a one unit change of $b_i$ results in $p_i$ units change of revenue. - -Thus, a dual variable can be interpreted as the **value** of one unit of resource $i$. - -This is why it is often called the **shadow price** of resource $i$. - -For feasible but not optimal primal and dual solutions $x$ and $p$, by weak duality, we have - -$$ -c_1 x_1 + c_2 x_2 + \dots + c_n x_n < b_1 p_1 + b_2 p_2 + \dots + b_m p_m. -$$ - -```{note} -Here, the expression is opposite to the statement above since primal problem is a minimization problem. -``` - -When a strict inequality holds, the solution is not optimal because it doesn't fully utilize all valuable resources. - -Evidently, - -* if a shadow price $p_i$ is larger than the market price for Resource $i$, the factory should buy more Resource $i$ and expand its scale to generate more revenue; - -* if a shadow price $p_i$ is less than the market price for Resource $i$, the factory should sell its Resource $i$. - -**Complementary slackness:** If there exists $i$ such that $a_i' x - b_i < 0$ for some $i$, then $p_i = 0$ by complementary slackness. $a_i' x - b_i < 0$ means that to achieve its optimal production, the factory doesn't require as much Resource $i$ as it has. -It is reasonable that the shadow price of Resource $i$ is 0: some of its resource $i$ is redundant. - -If there exists $j$ such that $A_j' p - c_j > 0$, then $x_j = 0$ by complementary slackness. $A_j' p - c_j > 0$ means that the value of all resources used when producing one unit of product $j$ is greater than its cost. - -This means that producing another product that can more efficiently utilize these resources is a better choice than producing product $j$ - -Since producing product $j$ is not optimal, $x_j$ should equal $0$. - -### Example 1: Production Problem - -This problem is one specific instance of the problem {eq}`linprog2`, whose economic meaning is interpreted above. - -Its dual problem is: - -$$ -\begin{aligned} -\min_{x_1,x_2} \ & 30 p_1 + 20 p_2 \\ -\mbox{subject to } \ & 2 p_1 + 4 p_2 \ge 3 \\ -& 5 p_1 + 2 p_2 \ge 4 \\ -& p_1, p_2 \ge 0 \\ -\end{aligned} -$$ - -We solve this dual problem by using the function *linprog*. - -Since parameters used here are defined before when solving the primal problem, we won't define them here. - -```{code-cell} ipython3 -# Solve the dual problem -res_ex1_dual = linprog(b_ex1, A_ub=-A_ex1.T, b_ub=-c_ex1) - -res_ex1_dual -``` - -The optimal value for the dual problem equals 27.5. - -This equals the optimal value of the primal problem, an illustration of strong duality. - -Shadow prices for materials and labor are 0.625 and 0.4375, respectively. - -### Example 2: Investment Problem - -The dual problem is: - -$$ -\begin{aligned} -\min_{p} \ & 100,000 p_1 - 20,000 p_4 - 20,000 p_5 - 20,000 p_6 + 50,000 p_7 \\ -\mbox{subject to } \ & p_1 + p_2 + p_3 \ge 1.30 \cdot 3 \\ -& p_1 - 1.06 p_2 + p_4 = 0 \\ -& p_2 - 1.06 p_3 + p_5 = 0 \\ -& p_3 + p_6 = 1.06 \\ -& p_2 + p_7 \ge 1.30 \\ -& p_i \ \text{unrestricted}, \quad i = 1,2,3 \\ -& p_i \le 0, \quad i = 4,5,6 \\ -& p_7 \ge 0 \\ -\end{aligned} -$$ - -We solve this dual problem by using the function *linprog*. - -```{code-cell} ipython3 -# Objective function parameters -c_ex2_dual = np.array([100000, 0, 0, -20000, -20000, -20000, 50000]) - -# Equality constraints -A_eq_ex2_dual = np.array([[1, -1.06, 0, 1, 0, 0, 0], - [0, 1, -1.06, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 1, 0]]) -b_eq_ex2_dual = np.array([0, 0, 1.06]) - -# Inequality constraints -A_ub_ex2_dual = - np.array([[1, 1, 1, 0, 0, 0, 0], - [0, 1, 0, 0, 0, 0, 1]]) -b_ub_ex2_dual = - np.array([1.30*3, 1.30]) - -# Bounds on decision variables -bounds_ex2_dual = [(None, None), - (None, None), - (None, None), - (None, 0), - (None, 0), - (None, 0), - ( 0, None)] - -# Solve the dual problem -res_ex2_dual = linprog(c_ex2_dual, A_eq=A_eq_ex2_dual, b_eq=b_eq_ex2_dual, - A_ub=A_ub_ex2_dual, b_ub=b_ub_ex2_dual, bounds=bounds_ex2_dual) - -res_ex2_dual -``` - -The optimal value for the dual problem is 141018.24, which equals the value of the primal problem. - -Now, let's interpret the dual variables. - -By strong duality and also our numerical results, we have that optimal value is: - -$$ -100,000 p_1 - 20,000 p_4 - 20,000 p_5 - 20,000 p_6 + 50,000 p_7. -$$ - -We know if $b_i$ changes one dollor, then the optimal payoff in the end of the third year will change $p_i$ dollars. - -For $i = 1$, this means if the initial capital changes by one dollar, then the optimal payoff in the end of the third year will change $p_1$ dollars. - -Thus, $p_1$ is the potential value of one more unit of initial capital, or the shadow price for initial capital. - -We can also interpret $p_1$ as the prospective value in the end of the third year coming from having one more dollar to invest at the beginning of the first year. - -If the mutual fund can raise money at a cost lower than $p_1 - 1$, then it should raise more money to increase its revenue. - -But if it bears a cost of funds higher than $p_1 - 1$, the mutual fund shouldn't do that. - -For $i = 4, 5, 6$, this means that if the amount of capital that the fund is permitted to borrow from the bank changes by one dollar, the optimal pay out at the end of the third year will change $p_i$ dollars. - -Thus, for $i = 4, 5, 6$, $|p_i|$ indicates the value of one dollar that the mutual fund can borrow from the bank at the beginning of the $i-3$-th year. - -$|p_i|$ is the shadow price for the loan amount. (We use absolute value here since $p_i \le 0$.) - -If the interest rate is lower than $|p_i|$, then the mutual fund should borrow to increase its optimal payoff; if the interest rate is higher, it is better to not do this. - -For $i = 7$, this means that if the amount of the corporate bond the mutual fund can buy changes one dollar, then the optimal payoff will change $p_7$ dollars at the end of the third year. Again, $p_7$ is the shadow price for the amount of the corporate bond the mutual fund can buy. - -As for numerical results - -1. $p_1 = 1.38$, which means one dollar of initial capital is worth $\$ 1.38$ at the end of the third year. - -2. $p_4 = p_5 = 0$, which means the loan amounts at the beginning of the first and second year are worth nothing. Recall that the optimal solution to the primal problem, $x_2, x_3 > 0$, which means at the beginning of the first and second year, the mutual fund has a postive bank account and borrows no capital from the bank. Thus, it is reasonable that the loan amounts at the beginning of the first and second year are valueless. This is what the complementary slackness conditions mean in this setting. - -3. $p_6 = -0.16$, which means one dollar of the loan amount at the beginning of the third year is worth $\$ 0.16$. Since $|p_6|$ is higher than the interest rate 6\%, the mutual fund should borrow as much as possible at the beginning of the third year. Recall that the optimal solution to the primal problem is $x_4 = -20,000$ which means the mutual fund borrows money from the bank as much as it can. - -4. $p_7 = 0.0015$, which means one dollar of the amount of the corporate bond that the mutual fund can buy is worth $\$ 0.0015$. diff --git a/lectures/short_path.md b/lectures/short_path.md deleted file mode 100644 index 13bd06d0e..000000000 --- a/lectures/short_path.md +++ /dev/null @@ -1,487 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -(short_path)= -```{raw} html -
- - QuantEcon - -
-``` - -# Shortest Paths - -```{index} single: Dynamic Programming; Shortest Paths -``` - -```{contents} Contents -:depth: 2 -``` - -## Overview - -The shortest path problem is a [classic problem](https://en.wikipedia.org/wiki/Shortest_path) in mathematics and computer science with applications in - -* Economics (sequential decision making, analysis of social networks, etc.) -* Operations research and transportation -* Robotics and artificial intelligence -* Telecommunication network design and routing -* etc., etc. - -Variations of the methods we discuss in this lecture are used millions of times every day, in applications such as - -* Google Maps -* routing packets on the internet - -For us, the shortest path problem also provides a nice introduction to the logic of **dynamic programming**. - -Dynamic programming is an extremely powerful optimization technique that we apply in many lectures on this site. - -The only scientific library we'll need in what follows is NumPy: - -```{code-cell} python3 -import numpy as np -``` - -## Outline of the Problem - -The shortest path problem is one of finding how to traverse a [graph](https://en.wikipedia.org/wiki/Graph_%28mathematics%29) from one specified node to another at minimum cost. - -Consider the following graph - -```{figure} /_static/lecture_specific/short_path/graph.png - -``` - -We wish to travel from node (vertex) A to node G at minimum cost - -* Arrows (edges) indicate the movements we can take. -* Numbers on edges indicate the cost of traveling that edge. - -(Graphs such as the one above are called weighted [directed graphs](https://en.wikipedia.org/wiki/Directed_graph).) - -Possible interpretations of the graph include - -* Minimum cost for supplier to reach a destination. -* Routing of packets on the internet (minimize time). -* Etc., etc. - -For this simple graph, a quick scan of the edges shows that the optimal paths are - -* A, C, F, G at cost 8 - -```{figure} /_static/lecture_specific/short_path/graph4.png - -``` - -* A, D, F, G at cost 8 - -```{figure} /_static/lecture_specific/short_path/graph3.png - -``` - -## Finding Least-Cost Paths - -For large graphs, we need a systematic solution. - -Let $J(v)$ denote the minimum cost-to-go from node $v$, understood as the total cost from $v$ if we take the best route. - -Suppose that we know $J(v)$ for each node $v$, as shown below for the graph from the preceding example - -```{figure} /_static/lecture_specific/short_path/graph2.png - -``` - -Note that $J(G) = 0$. - -The best path can now be found as follows - -1. Start at node $v = A$ -1. From current node $v$, move to any node that solves - -```{math} -:label: spprebell - -\min_{w \in F_v} \{ c(v, w) + J(w) \} -``` - -where - -* $F_v$ is the set of nodes that can be reached from $v$ in one step. -* $c(v, w)$ is the cost of traveling from $v$ to $w$. - -Hence, if we know the function $J$, then finding the best path is almost trivial. - -But how can we find the cost-to-go function $J$? - -Some thought will convince you that, for every node $v$, -the function $J$ satisfies - -```{math} -:label: spbell - -J(v) = \min_{w \in F_v} \{ c(v, w) + J(w) \} -``` - -This is known as the *Bellman equation*, after the mathematician Richard Bellman. - -The Bellman equation can be thought of as a restriction that $J$ must -satisfy. - -What we want to do now is use this restriction to compute $J$. - -## Solving for Minimum Cost-to-Go - -Let's look at an algorithm for computing $J$ and then think about how to -implement it. - -### The Algorithm - -The standard algorithm for finding $J$ is to start an initial guess and then iterate. - -This is a standard approach to solving nonlinear equations, often called -the method of **successive approximations**. - -Our initial guess will be - -```{math} -:label: spguess - -J_0(v) = 0 \text{ for all } v -``` - -Now - -1. Set $n = 0$ -1. Set $J_{n+1} (v) = \min_{w \in F_v} \{ c(v, w) + J_n(w) \}$ for all $v$ -1. If $J_{n+1}$ and $J_n$ are not equal then increment $n$, go to 2 - -This sequence converges to $J$. - -Although we omit the proof, we'll prove similar claims in our other lectures -on dynamic programming. - -### Implementation - -Having an algorithm is a good start, but we also need to think about how to -implement it on a computer. - -First, for the cost function $c$, we'll implement it as a matrix -$Q$, where a typical element is - -$$ -Q(v, w) -= -\begin{cases} - & c(v, w) \text{ if } w \in F_v \\ - & +\infty \text{ otherwise } -\end{cases} -$$ - -In this context $Q$ is usually called the **distance matrix**. - -We're also numbering the nodes now, with $A = 0$, so, for example - -$$ -Q(1, 2) -= -\text{ the cost of traveling from B to C } -$$ - -For example, for the simple graph above, we set - -```{code-cell} python3 -from numpy import inf - -Q = np.array([[inf, 1, 5, 3, inf, inf, inf], - [inf, inf, inf, 9, 6, inf, inf], - [inf, inf, inf, inf, inf, 2, inf], - [inf, inf, inf, inf, inf, 4, 8], - [inf, inf, inf, inf, inf, inf, 4], - [inf, inf, inf, inf, inf, inf, 1], - [inf, inf, inf, inf, inf, inf, 0]]) -``` - -Notice that the cost of staying still (on the principle diagonal) is set to - -* np.inf for non-destination nodes --- moving on is required. -* 0 for the destination node --- here is where we stop. - -For the sequence of approximations $\{J_n\}$ of the cost-to-go functions, we can use NumPy arrays. - -Let's try with this example and see how we go: - -```{code-cell} python3 -nodes = range(7) # Nodes = 0, 1, ..., 6 -J = np.zeros_like(nodes, dtype=int) # Initial guess -next_J = np.empty_like(nodes, dtype=int) # Stores updated guess - -max_iter = 500 -i = 0 - -while i < max_iter: - for v in nodes: - # minimize Q[v, w] + J[w] over all choices of w - lowest_cost = inf - for w in nodes: - cost = Q[v, w] + J[w] - if cost < lowest_cost: - lowest_cost = cost - next_J[v] = lowest_cost - if np.equal(next_J, J).all(): - break - else: - J[:] = next_J # Copy contents of next_J to J - i += 1 - -print("The cost-to-go function is", J) -``` - -This matches with the numbers we obtained by inspection above. - -But, importantly, we now have a methodology for tackling large graphs. - -## Exercises - - -```{exercise-start} -:label: short_path_ex1 -``` - -The text below describes a weighted directed graph. - -The line `node0, node1 0.04, node8 11.11, node14 72.21` means that from node0 we can go to - -* node1 at cost 0.04 -* node8 at cost 11.11 -* node14 at cost 72.21 - -No other nodes can be reached directly from node0. - -Other lines have a similar interpretation. - -Your task is to use the algorithm given above to find the optimal path and its cost. - -```{note} -You will be dealing with floating point numbers now, rather than -integers, so consider replacing `np.equal()` with `np.allclose()`. -``` - -```{code-cell} python3 -%%file graph.txt -node0, node1 0.04, node8 11.11, node14 72.21 -node1, node46 1247.25, node6 20.59, node13 64.94 -node2, node66 54.18, node31 166.80, node45 1561.45 -node3, node20 133.65, node6 2.06, node11 42.43 -node4, node75 3706.67, node5 0.73, node7 1.02 -node5, node45 1382.97, node7 3.33, node11 34.54 -node6, node31 63.17, node9 0.72, node10 13.10 -node7, node50 478.14, node9 3.15, node10 5.85 -node8, node69 577.91, node11 7.45, node12 3.18 -node9, node70 2454.28, node13 4.42, node20 16.53 -node10, node89 5352.79, node12 1.87, node16 25.16 -node11, node94 4961.32, node18 37.55, node20 65.08 -node12, node84 3914.62, node24 34.32, node28 170.04 -node13, node60 2135.95, node38 236.33, node40 475.33 -node14, node67 1878.96, node16 2.70, node24 38.65 -node15, node91 3597.11, node17 1.01, node18 2.57 -node16, node36 392.92, node19 3.49, node38 278.71 -node17, node76 783.29, node22 24.78, node23 26.45 -node18, node91 3363.17, node23 16.23, node28 55.84 -node19, node26 20.09, node20 0.24, node28 70.54 -node20, node98 3523.33, node24 9.81, node33 145.80 -node21, node56 626.04, node28 36.65, node31 27.06 -node22, node72 1447.22, node39 136.32, node40 124.22 -node23, node52 336.73, node26 2.66, node33 22.37 -node24, node66 875.19, node26 1.80, node28 14.25 -node25, node70 1343.63, node32 36.58, node35 45.55 -node26, node47 135.78, node27 0.01, node42 122.00 -node27, node65 480.55, node35 48.10, node43 246.24 -node28, node82 2538.18, node34 21.79, node36 15.52 -node29, node64 635.52, node32 4.22, node33 12.61 -node30, node98 2616.03, node33 5.61, node35 13.95 -node31, node98 3350.98, node36 20.44, node44 125.88 -node32, node97 2613.92, node34 3.33, node35 1.46 -node33, node81 1854.73, node41 3.23, node47 111.54 -node34, node73 1075.38, node42 51.52, node48 129.45 -node35, node52 17.57, node41 2.09, node50 78.81 -node36, node71 1171.60, node54 101.08, node57 260.46 -node37, node75 269.97, node38 0.36, node46 80.49 -node38, node93 2767.85, node40 1.79, node42 8.78 -node39, node50 39.88, node40 0.95, node41 1.34 -node40, node75 548.68, node47 28.57, node54 53.46 -node41, node53 18.23, node46 0.28, node54 162.24 -node42, node59 141.86, node47 10.08, node72 437.49 -node43, node98 2984.83, node54 95.06, node60 116.23 -node44, node91 807.39, node46 1.56, node47 2.14 -node45, node58 79.93, node47 3.68, node49 15.51 -node46, node52 22.68, node57 27.50, node67 65.48 -node47, node50 2.82, node56 49.31, node61 172.64 -node48, node99 2564.12, node59 34.52, node60 66.44 -node49, node78 53.79, node50 0.51, node56 10.89 -node50, node85 251.76, node53 1.38, node55 20.10 -node51, node98 2110.67, node59 23.67, node60 73.79 -node52, node94 1471.80, node64 102.41, node66 123.03 -node53, node72 22.85, node56 4.33, node67 88.35 -node54, node88 967.59, node59 24.30, node73 238.61 -node55, node84 86.09, node57 2.13, node64 60.80 -node56, node76 197.03, node57 0.02, node61 11.06 -node57, node86 701.09, node58 0.46, node60 7.01 -node58, node83 556.70, node64 29.85, node65 34.32 -node59, node90 820.66, node60 0.72, node71 0.67 -node60, node76 48.03, node65 4.76, node67 1.63 -node61, node98 1057.59, node63 0.95, node64 4.88 -node62, node91 132.23, node64 2.94, node76 38.43 -node63, node66 4.43, node72 70.08, node75 56.34 -node64, node80 47.73, node65 0.30, node76 11.98 -node65, node94 594.93, node66 0.64, node73 33.23 -node66, node98 395.63, node68 2.66, node73 37.53 -node67, node82 153.53, node68 0.09, node70 0.98 -node68, node94 232.10, node70 3.35, node71 1.66 -node69, node99 247.80, node70 0.06, node73 8.99 -node70, node76 27.18, node72 1.50, node73 8.37 -node71, node89 104.50, node74 8.86, node91 284.64 -node72, node76 15.32, node84 102.77, node92 133.06 -node73, node83 52.22, node76 1.40, node90 243.00 -node74, node81 1.07, node76 0.52, node78 8.08 -node75, node92 68.53, node76 0.81, node77 1.19 -node76, node85 13.18, node77 0.45, node78 2.36 -node77, node80 8.94, node78 0.98, node86 64.32 -node78, node98 355.90, node81 2.59 -node79, node81 0.09, node85 1.45, node91 22.35 -node80, node92 121.87, node88 28.78, node98 264.34 -node81, node94 99.78, node89 39.52, node92 99.89 -node82, node91 47.44, node88 28.05, node93 11.99 -node83, node94 114.95, node86 8.75, node88 5.78 -node84, node89 19.14, node94 30.41, node98 121.05 -node85, node97 94.51, node87 2.66, node89 4.90 -node86, node97 85.09 -node87, node88 0.21, node91 11.14, node92 21.23 -node88, node93 1.31, node91 6.83, node98 6.12 -node89, node97 36.97, node99 82.12 -node90, node96 23.53, node94 10.47, node99 50.99 -node91, node97 22.17 -node92, node96 10.83, node97 11.24, node99 34.68 -node93, node94 0.19, node97 6.71, node99 32.77 -node94, node98 5.91, node96 2.03 -node95, node98 6.17, node99 0.27 -node96, node98 3.32, node97 0.43, node99 5.87 -node97, node98 0.30 -node98, node99 0.33 -node99, -``` - -```{exercise-end} -``` - -```{solution-start} short_path_ex1 -:class: dropdown -``` - -First let's write a function that reads in the graph data above and builds a distance matrix. - -```{code-cell} python3 -num_nodes = 100 -destination_node = 99 - -def map_graph_to_distance_matrix(in_file): - - # First let's set of the distance matrix Q with inf everywhere - Q = np.full((num_nodes, num_nodes), np.inf) - - # Now we read in the data and modify Q - with open(in_file) as infile: - for line in infile: - elements = line.split(',') - node = elements.pop(0) - node = int(node[4:]) # convert node description to integer - if node != destination_node: - for element in elements: - destination, cost = element.split() - destination = int(destination[4:]) - Q[node, destination] = float(cost) - Q[destination_node, destination_node] = 0 - return Q -``` - -In addition, let's write - -1. a "Bellman operator" function that takes a distance matrix and current guess of J and returns an updated guess of J, and -1. a function that takes a distance matrix and returns a cost-to-go function. - -We'll use the algorithm described above. - -The minimization step is vectorized to make it faster. - -```{code-cell} python3 -def bellman(J, Q): - num_nodes = Q.shape[0] - next_J = np.empty_like(J) - for v in range(num_nodes): - next_J[v] = np.min(Q[v, :] + J) - return next_J - - -def compute_cost_to_go(Q): - num_nodes = Q.shape[0] - J = np.zeros(num_nodes) # Initial guess - max_iter = 500 - i = 0 - - while i < max_iter: - next_J = bellman(J, Q) - if np.allclose(next_J, J): - break - else: - J[:] = next_J # Copy contents of next_J to J - i += 1 - - return(J) -``` - -We used np.allclose() rather than testing exact equality because we are -dealing with floating point numbers now. - -Finally, here's a function that uses the cost-to-go function to obtain the -optimal path (and its cost). - -```{code-cell} python3 -def print_best_path(J, Q): - sum_costs = 0 - current_node = 0 - while current_node != destination_node: - print(current_node) - # Move to the next node and increment costs - next_node = np.argmin(Q[current_node, :] + J) - sum_costs += Q[current_node, next_node] - current_node = next_node - - print(destination_node) - print('Cost: ', sum_costs) -``` - -Okay, now we have the necessary functions, let's call them to do the job we were assigned. - -```{code-cell} python3 -Q = map_graph_to_distance_matrix('graph.txt') -J = compute_cost_to_go(Q) -print_best_path(J, Q) -``` - -The total cost of the path should agree with $J[0]$ so let's check this. - -```{code-cell} python3 -J[0] -``` - -```{solution-end} -``` \ No newline at end of file From b315f6b65707b907f8e4eff3f881deeee27521e7 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 29 Jan 2024 18:41:31 +1100 Subject: [PATCH 2/5] update toc and tags --- lectures/_toc.yml | 3 --- lectures/cake_eating_problem.md | 2 +- lectures/lqcontrol.md | 4 ++-- lectures/opt_transport.md | 4 ++-- lectures/scalar_dynam.md | 2 +- 5 files changed, 6 insertions(+), 9 deletions(-) diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 588614dc1..451de3804 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -4,7 +4,6 @@ parts: - caption: Tools and Techniques numbered: true chapters: - - file: geom_series - file: sir_model - file: linear_algebra - file: qr_decomp @@ -31,7 +30,6 @@ parts: - caption: Linear Programming numbered: true chapters: - - file: lp_intro - file: opt_transport - file: von_neumann_model @@ -47,7 +45,6 @@ parts: - file: wealth_dynamics - file: kalman - file: kalman_2 - - file: short_path - caption: Search numbered: true chapters: diff --git a/lectures/cake_eating_problem.md b/lectures/cake_eating_problem.md index 7ba80b36d..a4627f706 100644 --- a/lectures/cake_eating_problem.md +++ b/lectures/cake_eating_problem.md @@ -32,7 +32,7 @@ The main tool we will use to solve the cake eating problem is dynamic programmin Readers might find it helpful to review the following lectures before reading this one: -* The {doc}`shortest paths lecture ` +* The {doc}`shortest paths lecture ` * The {doc}`basic McCall model ` * The {doc}`McCall model with separation ` * The {doc}`McCall model with separation and a continuous wage distribution ` diff --git a/lectures/lqcontrol.md b/lectures/lqcontrol.md index ce3336a68..6169ad61d 100644 --- a/lectures/lqcontrol.md +++ b/lectures/lqcontrol.md @@ -57,7 +57,7 @@ In reading what follows, it will be useful to have some familiarity with * matrix manipulations * vectors of random variables -* dynamic programming and the Bellman equation (see for example {doc}`this lecture ` and {doc}`this lecture `) +* dynamic programming and the Bellman equation (see for example {doc}`this lecture ` and {doc}`this lecture `) For additional reading on LQ control, see, for example, @@ -369,7 +369,7 @@ What's special about the LQ case is that -- as we shall soon see --- the optima ### Solution To solve the finite horizon LQ problem we can use a dynamic programming -strategy based on backward induction that is conceptually similar to the approach adopted in {doc}`this lecture `. +strategy based on backward induction that is conceptually similar to the approach adopted in {doc}`this lecture `. For reasons that will soon become clear, we first introduce the notation $J_T(x) = x' R_f x$. diff --git a/lectures/opt_transport.md b/lectures/opt_transport.md index 5770df5e6..e6d1ab348 100644 --- a/lectures/opt_transport.md +++ b/lectures/opt_transport.md @@ -20,7 +20,7 @@ because of its many applications and because of its important role in the histor economic theory. In this lecture, we describe the problem, tell how -{doc}`linear programming ` is a +{doc}`linear programming ` is a key tool for solving it, and then provide some examples. We will provide other applications in followup lectures. @@ -631,7 +631,7 @@ linprog(-b[:-1], A_ub=A[:-1].T, b_ub=C_vec, ### Interpretation of dual problem By **strong duality** (please see this lecture -{doc}`Linear Programming `), we know that: +{doc}`Linear Programming `), we know that: $$ \sum_{i=1}^m \sum_{j=1}^n c_{ij} x_{ij} = \sum_{i=1}^m p_i u_i + \sum_{j=1}^n q_j v_j diff --git a/lectures/scalar_dynam.md b/lectures/scalar_dynam.md index 7a78b83e9..f09f1f9c4 100644 --- a/lectures/scalar_dynam.md +++ b/lectures/scalar_dynam.md @@ -100,7 +100,7 @@ a x_0 + b, \quad a^2 x_0 + a b + b, \quad \text{etc.} ``` -Continuing in this way, and using our knowledge of {doc}`geometric series `, we find that, for any $t \geq 0$, +Continuing in this way, and using our knowledge of {doc}`geometric series `, we find that, for any $t \geq 0$, ```{math} :label: sdslinmod From 7bcd2d716fccfb49a742ac4289bd5de844057a2f Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 29 Jan 2024 19:32:37 +1100 Subject: [PATCH 3/5] remove schelling --- lectures/_toc.yml | 1 - lectures/schelling.md | 279 ------------------------------------------ 2 files changed, 280 deletions(-) delete mode 100644 lectures/schelling.md diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 451de3804..1e1723fc0 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -100,7 +100,6 @@ parts: - caption: Multiple Agent Models numbered: true chapters: - - file: schelling - file: lake_model - file: rational_expectations - file: re_with_feedback diff --git a/lectures/schelling.md b/lectures/schelling.md deleted file mode 100644 index 6fdf06253..000000000 --- a/lectures/schelling.md +++ /dev/null @@ -1,279 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -(schelling)= -```{raw} html - -``` - -# Schelling's Segregation Model - -```{index} single: Schelling Segregation Model -``` - -```{index} single: Models; Schelling's Segregation Model -``` - -```{contents} Contents -:depth: 2 -``` - -## Outline - -In 1969, Thomas C. Schelling developed a simple but striking model of racial segregation {cite}`Schelling1969`. - -His model studies the dynamics of racially mixed neighborhoods. - -Like much of Schelling's work, the model shows how local interactions can lead to surprising aggregate structure. - -In particular, it shows that relatively mild preference for neighbors of similar race can lead in aggregate to the collapse of mixed neighborhoods, and high levels of segregation. - -In recognition of this and other research, Schelling was awarded the 2005 Nobel Prize in Economic Sciences (joint with Robert Aumann). - -In this lecture, we (in fact you) will build and run a version of Schelling's model. - -Let's start with some imports: - -```{code-cell} ipython -import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size -from random import uniform, seed -from math import sqrt -``` - -## The Model - -We will cover a variation of Schelling's model that is easy to program and captures the main idea. - -### Set-Up - -Suppose we have two types of people: orange people and green people. - -For the purpose of this lecture, we will assume there are 250 of each type. - -These agents all live on a single unit square. - -The location of an agent is just a point $(x, y)$, where $0 < x, y < 1$. - -### Preferences - -We will say that an agent is *happy* if half or more of her 10 nearest neighbors are of the same type. - -Here 'nearest' is in terms of [Euclidean distance](https://en.wikipedia.org/wiki/Euclidean_distance). - -An agent who is not happy is called *unhappy*. - -An important point here is that agents are not averse to living in mixed areas. - -They are perfectly happy if half their neighbors are of the other color. - -### Behavior - -Initially, agents are mixed together (integrated). - -In particular, the initial location of each agent is an independent draw from a bivariate uniform distribution on $S = (0, 1)^2$. - -Now, cycling through the set of all agents, each agent is now given the chance to stay or move. - -We assume that each agent will stay put if they are happy and move if unhappy. - -The algorithm for moving is as follows - -1. Draw a random location in $S$ -1. If happy at new location, move there -1. Else, go to step 1 - -In this way, we cycle continuously through the agents, moving as required. - -We continue to cycle until no one wishes to move. - -## Results - -Let's have a look at the results we got when we coded and ran this model. - -As discussed above, agents are initially mixed randomly together. - -```{figure} /_static/lecture_specific/schelling/schelling_fig1.png - -``` - -But after several cycles, they become segregated into distinct regions. - -```{figure} /_static/lecture_specific/schelling/schelling_fig2.png - -``` - -```{figure} /_static/lecture_specific/schelling/schelling_fig3.png - -``` - -```{figure} /_static/lecture_specific/schelling/schelling_fig4.png - -``` - -In this instance, the program terminated after 4 cycles through the set of -agents, indicating that all agents had reached a state of happiness. - -What is striking about the pictures is how rapidly racial integration breaks down. - -This is despite the fact that people in the model don't actually mind living mixed with the other type. - -Even with these preferences, the outcome is a high degree of segregation. - -## Exercises - -```{exercise-start} -:label: schelling_ex1 -``` - -Implement and run this simulation for yourself. - -Consider the following structure for your program. - -Agents can be modeled as [objects](https://python-programming.quantecon.org/python_oop.html). - -Here's an indication of how they might look - -```{code-block} none -* Data: - - * type (green or orange) - * location - -* Methods: - - * determine whether happy or not given locations of other agents - - * If not happy, move - - * find a new location where happy -``` - -And here's some pseudocode for the main loop - -```{code-block} none -while agents are still moving - for agent in agents - give agent the opportunity to move -``` - -Use 250 agents of each type. - -```{exercise-end} -``` - -```{solution-start} schelling_ex1 -:class: dropdown -``` - -Here's one solution that does the job we want. - -If you feel like a further exercise, you can probably speed up some of the computations and -then increase the number of agents. - -```{code-cell} python3 -seed(10) # For reproducible random numbers - -class Agent: - - def __init__(self, type): - self.type = type - self.draw_location() - - def draw_location(self): - self.location = uniform(0, 1), uniform(0, 1) - - def get_distance(self, other): - "Computes the euclidean distance between self and other agent." - a = (self.location[0] - other.location[0])**2 - b = (self.location[1] - other.location[1])**2 - return sqrt(a + b) - - def happy(self, agents): - "True if sufficient number of nearest neighbors are of the same type." - distances = [] - # distances is a list of pairs (d, agent), where d is distance from - # agent to self - for agent in agents: - if self != agent: - distance = self.get_distance(agent) - distances.append((distance, agent)) - # == Sort from smallest to largest, according to distance == # - distances.sort() - # == Extract the neighboring agents == # - neighbors = [agent for d, agent in distances[:num_neighbors]] - # == Count how many neighbors have the same type as self == # - num_same_type = sum(self.type == agent.type for agent in neighbors) - return num_same_type >= require_same_type - - def update(self, agents): - "If not happy, then randomly choose new locations until happy." - while not self.happy(agents): - self.draw_location() - - -def plot_distribution(agents, cycle_num): - "Plot the distribution of agents after cycle_num rounds of the loop." - x_values_0, y_values_0 = [], [] - x_values_1, y_values_1 = [], [] - # == Obtain locations of each type == # - for agent in agents: - x, y = agent.location - if agent.type == 0: - x_values_0.append(x) - y_values_0.append(y) - else: - x_values_1.append(x) - y_values_1.append(y) - fig, ax = plt.subplots(figsize=(8, 8)) - plot_args = {'markersize': 8, 'alpha': 0.6} - ax.set_facecolor('azure') - ax.plot(x_values_0, y_values_0, 'o', markerfacecolor='orange', **plot_args) - ax.plot(x_values_1, y_values_1, 'o', markerfacecolor='green', **plot_args) - ax.set_title(f'Cycle {cycle_num-1}') - plt.show() - -# == Main == # - -num_of_type_0 = 250 -num_of_type_1 = 250 -num_neighbors = 10 # Number of agents regarded as neighbors -require_same_type = 5 # Want at least this many neighbors to be same type - -# == Create a list of agents == # -agents = [Agent(0) for i in range(num_of_type_0)] -agents.extend(Agent(1) for i in range(num_of_type_1)) - - -count = 1 -# == Loop until none wishes to move == # -while True: - print('Entering loop ', count) - plot_distribution(agents, count) - count += 1 - no_one_moved = True - for agent in agents: - old_location = agent.location - agent.update(agents) - if agent.location != old_location: - no_one_moved = False - if no_one_moved: - break - -print('Converged, terminating.') -``` - -```{solution-end} -``` From d42e8d2a2f9b4bb710302d3fb2f99dd9fe64eab5 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 29 Jan 2024 19:43:12 +1100 Subject: [PATCH 4/5] add schelling redirect --- lectures/_config.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/lectures/_config.yml b/lectures/_config.yml index baf81c18f..ee7ff19b7 100644 --- a/lectures/_config.yml +++ b/lectures/_config.yml @@ -107,6 +107,7 @@ sphinx: geom_series: https://intro.quantecon.org/geom_series.html lp_intro: https://intro.quantecon.org/lp_intro.html short_path: https://intro.quantecon.org/short_path.html + schelling: https://intro.quantecon.org/schelling.html tojupyter_static_file_path: ["source/_static", "_static"] tojupyter_target_html: true tojupyter_urlpath: "https://python.quantecon.org/" From aba6fa53a70e6beb5d3c92b58c0adac497e752f3 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Wed, 31 Jan 2024 22:24:36 +1100 Subject: [PATCH 5/5] migrate scalar_dynam and complex_and_trig --- lectures/_config.yml | 2 + lectures/_toc.yml | 3 - lectures/complex_and_trig.md | 528 ----------------------------------- lectures/heavy_tails.md | 28 -- lectures/scalar_dynam.md | 518 ---------------------------------- 5 files changed, 2 insertions(+), 1077 deletions(-) delete mode 100644 lectures/complex_and_trig.md delete mode 100644 lectures/heavy_tails.md delete mode 100644 lectures/scalar_dynam.md diff --git a/lectures/_config.yml b/lectures/_config.yml index ee7ff19b7..da95b8551 100644 --- a/lectures/_config.yml +++ b/lectures/_config.yml @@ -108,6 +108,8 @@ sphinx: lp_intro: https://intro.quantecon.org/lp_intro.html short_path: https://intro.quantecon.org/short_path.html schelling: https://intro.quantecon.org/schelling.html + scalar_dynam: https://intro.quantecon.org/scalar_dynam.html + complex_and_trig: https://intro.quantecon.org/complex_and_trig.html tojupyter_static_file_path: ["source/_static", "_static"] tojupyter_target_html: true tojupyter_urlpath: "https://python.quantecon.org/" diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 1e1723fc0..a4f41e6f7 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -7,7 +7,6 @@ parts: - file: sir_model - file: linear_algebra - file: qr_decomp - - file: complex_and_trig - file: eig_circulant - file: svd_intro - file: var_dmd @@ -21,7 +20,6 @@ parts: - file: prob_meaning - file: multi_hyper - file: multivariate_normal - - file: heavy_tails - file: hoist_failure - file: back_prop - file: rand_resp @@ -36,7 +34,6 @@ parts: - caption: Introduction to Dynamics numbered: true chapters: - - file: scalar_dynam - file: finite_markov - file: inventory_dynamics - file: linear_models diff --git a/lectures/complex_and_trig.md b/lectures/complex_and_trig.md deleted file mode 100644 index 92ebb783b..000000000 --- a/lectures/complex_and_trig.md +++ /dev/null @@ -1,528 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -(complex_and_trig)= -```{raw} html - -``` - -```{index} single: python -``` - -# Complex Numbers and Trigonometry - -```{contents} Contents -:depth: 2 -``` - -## Overview - -This lecture introduces some elementary mathematics and trigonometry. - -Useful and interesting in its own right, these concepts reap substantial rewards when studying dynamics generated -by linear difference equations or linear differential equations. - -For example, these tools are keys to understanding outcomes attained by Paul -Samuelson (1939) {cite}`Samuelson1939` in his classic paper on interactions -between the investment accelerator and the Keynesian consumption function, our -topic in the lecture {doc}`Samuelson Multiplier Accelerator `. - -In addition to providing foundations for Samuelson's work and extensions of -it, this lecture can be read as a stand-alone quick reminder of key results -from elementary high school trigonometry. - -So let's dive in. - -### Complex Numbers - -A complex number has a **real part** $x$ and a purely **imaginary part** $y$. - -The Euclidean, polar, and trigonometric forms of a complex number $z$ are: - -$$ -z = x + iy = re^{i\theta} = r(\cos{\theta} + i \sin{\theta}) -$$ - -The second equality above is known as **Euler's formula** - -- [Euler](https://en.wikipedia.org/wiki/Leonhard_Euler) contributed many other formulas too! - -The complex conjugate $\bar z$ of $z$ is defined as - -$$ -\bar z = x - iy = r e^{-i \theta} = r (\cos{\theta} - i \sin{\theta} ) -$$ - -The value $x$ is the **real** part of $z$ and $y$ is the -**imaginary** part of $z$. - -The symbol $| z |$ = $\sqrt{\bar{z}\cdot z} = r$ represents the **modulus** of $z$. - -The value $r$ is the Euclidean distance of vector $(x,y)$ from the -origin: - -$$ -r = |z| = \sqrt{x^2 + y^2} -$$ - -The value $\theta$ is the angle of $(x,y)$ with respect to the real axis. - -Evidently, the tangent of $\theta$ is $\left(\frac{y}{x}\right)$. - -Therefore, - -$$ -\theta = \tan^{-1} \Big( \frac{y}{x} \Big) -$$ - -Three elementary trigonometric functions are - -$$ -\cos{\theta} = \frac{x}{r} = \frac{e^{i\theta} + e^{-i\theta}}{2} , \quad -\sin{\theta} = \frac{y}{r} = \frac{e^{i\theta} - e^{-i\theta}}{2i} , \quad -\tan{\theta} = \frac{y}{x} -$$ - -We'll need the following imports: - -```{code-cell} ipython -import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size -import numpy as np -from sympy import (Symbol, symbols, Eq, nsolve, sqrt, cos, sin, simplify, - init_printing, integrate) -``` - -### An Example - -Consider the complex number $z = 1 + \sqrt{3} i$. - -For $z = 1 + \sqrt{3} i$, $x = 1$, $y = \sqrt{3}$. - -It follows that $r = 2$ and -$\theta = \tan^{-1}(\sqrt{3}) = \frac{\pi}{3} = 60^o$. - -Let's use Python to plot the trigonometric form of the complex number -$z = 1 + \sqrt{3} i$. - -```{code-cell} python3 -# Abbreviate useful values and functions -π = np.pi - - -# Set parameters -r = 2 -θ = π/3 -x = r * np.cos(θ) -x_range = np.linspace(0, x, 1000) -θ_range = np.linspace(0, θ, 1000) - -# Plot -fig = plt.figure(figsize=(8, 8)) -ax = plt.subplot(111, projection='polar') - -ax.plot((0, θ), (0, r), marker='o', color='b') # Plot r -ax.plot(np.zeros(x_range.shape), x_range, color='b') # Plot x -ax.plot(θ_range, x / np.cos(θ_range), color='b') # Plot y -ax.plot(θ_range, np.full(θ_range.shape, 0.1), color='r') # Plot θ - -ax.margins(0) # Let the plot starts at origin - -ax.set_title("Trigonometry of complex numbers", va='bottom', - fontsize='x-large') - -ax.set_rmax(2) -ax.set_rticks((0.5, 1, 1.5, 2)) # Less radial ticks -ax.set_rlabel_position(-88.5) # Get radial labels away from plotted line - -ax.text(θ, r+0.01 , r'$z = x + iy = 1 + \sqrt{3}\, i$') # Label z -ax.text(θ+0.2, 1 , '$r = 2$') # Label r -ax.text(0-0.2, 0.5, '$x = 1$') # Label x -ax.text(0.5, 1.2, r'$y = \sqrt{3}$') # Label y -ax.text(0.25, 0.15, r'$\theta = 60^o$') # Label θ - -ax.grid(True) -plt.show() -``` - -## De Moivre's Theorem - -de Moivre's theorem states that: - -$$ -(r(\cos{\theta} + i \sin{\theta}))^n = -r^n e^{in\theta} = -r^n(\cos{n\theta} + i \sin{n\theta}) -$$ - -To prove de Moivre's theorem, note that - -$$ -(r(\cos{\theta} + i \sin{\theta}))^n = \big( re^{i\theta} \big)^n -$$ - -and compute. - -## Applications of de Moivre's Theorem - -### Example 1 - -We can use de Moivre's theorem to show that -$r = \sqrt{x^2 + y^2}$. - -We have - -$$ -\begin{aligned} -1 &= e^{i\theta} e^{-i\theta} \\ -&= (\cos{\theta} + i \sin{\theta})(\cos{(\text{-}\theta)} + i \sin{(\text{-}\theta)}) \\ -&= (\cos{\theta} + i \sin{\theta})(\cos{\theta} - i \sin{\theta}) \\ -&= \cos^2{\theta} + \sin^2{\theta} \\ -&= \frac{x^2}{r^2} + \frac{y^2}{r^2} -\end{aligned} -$$ - -and thus - -$$ -x^2 + y^2 = r^2 -$$ - -We recognize this as a theorem of **Pythagoras**. - -### Example 2 - -Let $z = re^{i\theta}$ and $\bar{z} = re^{-i\theta}$ so that $\bar{z}$ is the **complex conjugate** of $z$. - -$(z, \bar z)$ form a **complex conjugate pair** of complex numbers. - -Let $a = pe^{i\omega}$ and $\bar{a} = pe^{-i\omega}$ be -another complex conjugate pair. - -For each element of a sequence of integers $n = 0, 1, 2, \ldots, $. - -To do so, we can apply de Moivre's formula. - -Thus, - -$$ -\begin{aligned} -x_n &= az^n + \bar{a}\bar{z}^n \\ -&= p e^{i\omega} (re^{i\theta})^n + p e^{-i\omega} (re^{-i\theta})^n \\ -&= pr^n e^{i (\omega + n\theta)} + pr^n e^{-i (\omega + n\theta)} \\ -&= pr^n [\cos{(\omega + n\theta)} + i \sin{(\omega + n\theta)} + - \cos{(\omega + n\theta)} - i \sin{(\omega + n\theta)}] \\ -&= 2 pr^n \cos{(\omega + n\theta)} -\end{aligned} -$$ - -### Example 3 - -This example provides machinery that is at the heard of Samuelson's analysis of his multiplier-accelerator model {cite}`Samuelson1939`. - -Thus, consider a **second-order linear difference equation** - -$$ -x_{n+2} = c_1 x_{n+1} + c_2 x_n -$$ - -whose **characteristic polynomial** is - -$$ -z^2 - c_1 z - c_2 = 0 -$$ - -or - -$$ -(z^2 - c_1 z - c_2 ) = (z - z_1)(z- z_2) = 0 -$$ - -has roots $z_1, z_1$. - -A **solution** is a sequence $\{x_n\}_{n=0}^\infty$ that satisfies -the difference equation. - -Under the following circumstances, we can apply our example 2 formula to -solve the difference equation - -- the roots $z_1, z_2$ of the characteristic polynomial of the - difference equation form a complex conjugate pair -- the values $x_0, x_1$ are given initial conditions - -To solve the difference equation, recall from example 2 that - -$$ -x_n = 2 pr^n \cos{(\omega + n\theta)} -$$ - -where $\omega, p$ are coefficients to be determined from -information encoded in the initial conditions $x_1, x_0$. - -Since -$x_0 = 2 p \cos{\omega}$ and $x_1 = 2 pr \cos{(\omega + \theta)}$ -the ratio of $x_1$ to $x_0$ is - -$$ -\frac{x_1}{x_0} = \frac{r \cos{(\omega + \theta)}}{\cos{\omega}} -$$ - -We can solve this equation for $\omega$ then solve for $p$ using $x_0 = 2 pr^0 \cos{(\omega + n\theta)}$. - -With the `sympy` package in Python, we are able to solve and plot the -dynamics of $x_n$ given different values of $n$. - -In this example, we set the initial values: - $r = 0.9$ - -$\theta = \frac{1}{4}\pi$ - $x_0 = 4$ - -$x_1 = r \cdot 2\sqrt{2} = 1.8 \sqrt{2}$. - -We first numerically solve for $\omega$ and $p$ using -`nsolve` in the `sympy` package based on the above initial -condition: - -```{code-cell} python3 -# Set parameters -r = 0.9 -θ = π/4 -x0 = 4 -x1 = 2 * r * sqrt(2) - -# Define symbols to be calculated -ω, p = symbols('ω p', real=True) - -# Solve for ω -## Note: we choose the solution near 0 -eq1 = Eq(x1/x0 - r * cos(ω+θ) / cos(ω), 0) -ω = nsolve(eq1, ω, 0) -ω = float(ω) -print(f'ω = {ω:1.3f}') - -# Solve for p -eq2 = Eq(x0 - 2 * p * cos(ω), 0) -p = nsolve(eq2, p, 0) -p = float(p) -print(f'p = {p:1.3f}') -``` - -Using the code above, we compute that -$\omega = 0$ and $p = 2$. - -Then we plug in the values we solve for $\omega$ and $p$ -and plot the dynamic. - -```{code-cell} python3 -# Define range of n -max_n = 30 -n = np.arange(0, max_n+1, 0.01) - -# Define x_n -x = lambda n: 2 * p * r**n * np.cos(ω + n * θ) - -# Plot -fig, ax = plt.subplots(figsize=(12, 8)) - -ax.plot(n, x(n)) -ax.set(xlim=(0, max_n), ylim=(-5, 5), xlabel='$n$', ylabel='$x_n$') - -# Set x-axis in the middle of the plot -ax.spines['bottom'].set_position('center') -ax.spines['right'].set_color('none') -ax.spines['top'].set_color('none') -ax.xaxis.set_ticks_position('bottom') -ax.yaxis.set_ticks_position('left') - -ticklab = ax.xaxis.get_ticklabels()[0] # Set x-label position -trans = ticklab.get_transform() -ax.xaxis.set_label_coords(31, 0, transform=trans) - -ticklab = ax.yaxis.get_ticklabels()[0] # Set y-label position -trans = ticklab.get_transform() -ax.yaxis.set_label_coords(0, 5, transform=trans) - -ax.grid() -plt.show() -``` - -### Trigonometric Identities - -We can obtain a complete suite of trigonometric identities by -appropriately manipulating polar forms of complex numbers. - -We'll get many of them by deducing implications of the equality - -$$ -e^{i(\omega + \theta)} = e^{i\omega} e^{i\theta} -$$ - -For example, we'll calculate identities for - -$\cos{(\omega + \theta)}$ and $\sin{(\omega + \theta)}$. - -Using the sine and cosine formulas presented at the beginning of this -lecture, we have: - -$$ -\begin{aligned} -\cos{(\omega + \theta)} = \frac{e^{i(\omega + \theta)} + e^{-i(\omega + \theta)}}{2} \\ -\sin{(\omega + \theta)} = \frac{e^{i(\omega + \theta)} - e^{-i(\omega + \theta)}}{2i} -\end{aligned} -$$ - -We can also obtain the trigonometric identities as follows: - -$$ -\begin{aligned} -\cos{(\omega + \theta)} + i \sin{(\omega + \theta)} -&= e^{i(\omega + \theta)} \\ -&= e^{i\omega} e^{i\theta} \\ -&= (\cos{\omega} + i \sin{\omega})(\cos{\theta} + i \sin{\theta}) \\ -&= (\cos{\omega}\cos{\theta} - \sin{\omega}\sin{\theta}) + -i (\cos{\omega}\sin{\theta} + \sin{\omega}\cos{\theta}) -\end{aligned} -$$ - -Since both real and imaginary parts of the above formula should be -equal, we get: - -$$ -\begin{aligned} -\cos{(\omega + \theta)} = \cos{\omega}\cos{\theta} - \sin{\omega}\sin{\theta} \\ -\sin{(\omega + \theta)} = \cos{\omega}\sin{\theta} + \sin{\omega}\cos{\theta} -\end{aligned} -$$ - -The equations above are also known as the **angle sum identities**. We -can verify the equations using the `simplify` function in the -`sympy` package: - -```{code-cell} python3 -# Define symbols -ω, θ = symbols('ω θ', real=True) - -# Verify -print("cos(ω)cos(θ) - sin(ω)sin(θ) =", - simplify(cos(ω)*cos(θ) - sin(ω) * sin(θ))) -print("cos(ω)sin(θ) + sin(ω)cos(θ) =", - simplify(cos(ω)*sin(θ) + sin(ω) * cos(θ))) -``` - -### Trigonometric Integrals - -We can also compute the trigonometric integrals using polar forms of -complex numbers. - -For example, we want to solve the following integral: - -$$ -\int_{-\pi}^{\pi} \cos(\omega) \sin(\omega) \, d\omega -$$ - -Using Euler's formula, we have: - -$$ -\begin{aligned} -\int \cos(\omega) \sin(\omega) \, d\omega -&= -\int -\frac{(e^{i\omega} + e^{-i\omega})}{2} -\frac{(e^{i\omega} - e^{-i\omega})}{2i} -\, d\omega \\ -&= -\frac{1}{4i} -\int -e^{2i\omega} - e^{-2i\omega} -\, d\omega \\ -&= -\frac{1}{4i} -\bigg( \frac{-i}{2} e^{2i\omega} - \frac{i}{2} e^{-2i\omega} + C_1 \bigg) \\ -&= --\frac{1}{8} -\bigg[ \bigg(e^{i\omega}\bigg)^2 + \bigg(e^{-i\omega}\bigg)^2 - 2 \bigg] + C_2 \\ -&= --\frac{1}{8} (e^{i\omega} - e^{-i\omega})^2 + C_2 \\ -&= -\frac{1}{2} \bigg( \frac{e^{i\omega} - e^{-i\omega}}{2i} \bigg)^2 + C_2 \\ -&= \frac{1}{2} \sin^2(\omega) + C_2 -\end{aligned} -$$ - -and thus: - -$$ -\int_{-\pi}^{\pi} \cos(\omega) \sin(\omega) \, d\omega = -\frac{1}{2}\sin^2(\pi) - \frac{1}{2}\sin^2(-\pi) = 0 -$$ - -We can verify the analytical as well as numerical results using -`integrate` in the `sympy` package: - -```{code-cell} python3 -# Set initial printing -init_printing() - -ω = Symbol('ω') -print('The analytical solution for integral of cos(ω)sin(ω) is:') -integrate(cos(ω) * sin(ω), ω) -``` - -```{code-cell} python3 -print('The numerical solution for the integral of cos(ω)sin(ω) \ -from -π to π is:') -integrate(cos(ω) * sin(ω), (ω, -π, π)) -``` - -### Exercises - -```{exercise} -:label: complex_ex1 - -We invite the reader to verify analytically and with the `sympy` package the following two equalities: - -$$ -\int_{-\pi}^{\pi} \cos (\omega)^2 \, d\omega = \pi -$$ - -$$ -\int_{-\pi}^{\pi} \sin (\omega)^2 \, d\omega = \pi -$$ -``` - -```{solution-start} complex_ex1 -:class: dropdown -``` - -Let's import symbolic $\pi$ from `sympy` - -```{code-cell} ipython3 -# Import symbolic π from sympy -from sympy import pi -``` - -```{code-cell} ipython3 -print('The analytical solution for the integral of cos(ω)**2 \ -from -π to π is:') - -integrate(cos(ω)**2, (ω, -pi, pi)) -``` - -```{code-cell} ipython3 -print('The analytical solution for the integral of sin(ω)**2 \ -from -π to π is:') - -integrate(sin(ω)**2, (ω, -pi, pi)) -``` - -```{solution-end} -``` diff --git a/lectures/heavy_tails.md b/lectures/heavy_tails.md deleted file mode 100644 index d42c1f40a..000000000 --- a/lectures/heavy_tails.md +++ /dev/null @@ -1,28 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -```{raw} html - -``` - -# Heavy-Tailed Distributions - -```{admonition} New website -:class: warning - -We have moved this lecture to a new lecture series. - -See [An Undergraduate Lecture Series for the Foundations of Computational Economics](https://intro.quantecon.org) -``` \ No newline at end of file diff --git a/lectures/scalar_dynam.md b/lectures/scalar_dynam.md deleted file mode 100644 index f09f1f9c4..000000000 --- a/lectures/scalar_dynam.md +++ /dev/null @@ -1,518 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -```{raw} html - -``` - -# {index}`Dynamics in One Dimension ` - -```{contents} Contents -:depth: 2 -``` - -## Overview - -In this lecture we give a quick introduction to discrete time dynamics in one -dimension. - -In one-dimensional models, the state of the system is described by a single variable. - -Although most interesting dynamic models have two or more state variables, the -one-dimensional setting is a good place to learn the foundations of dynamics and build -intuition. - -Let's start with some standard imports: - -```{code-cell} ipython -import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size -import numpy as np -``` - -## Some Definitions - -This section sets out the objects of interest and the kinds of properties we study. - -### Difference Equations - -A **time homogeneous first order difference equation** is an equation of the -form - -```{math} -:label: sdsod - -x_{t+1} = g(x_t) -``` - -where $g$ is a function from some subset $S$ of $\mathbb R$ to itself. - -Here $S$ is called the **state space** and $x$ is called the **state variable**. - -In the definition, - -* time homogeneity means that $g$ is the same at each time $t$ -* first order means dependence on only one lag (i.e., earlier states such as $x_{t-1}$ do not enter into {eq}`sdsod`). - -If $x_0 \in S$ is given, then {eq}`sdsod` recursively defines the sequence - -```{math} -:label: sdstraj - -x_0, \quad -x_1 = g(x_0), \quad -x_2 = g(x_1) = g(g(x_0)), \quad \text{etc.} -``` - -This sequence is called the **trajectory** of $x_0$ under $g$. - -If we define $g^n$ to be $n$ compositions of $g$ with itself, then we can write the trajectory more simply as $x_t = g^t(x_0)$ for $t \geq 0$. - -### Example: A Linear Model - -One simple example is the **linear difference equation** - -$$ -x_{t+1} = a x_t + b, \qquad S = \mathbb R -$$ - -where $a, b$ are fixed constants. - -In this case, given $x_0$, the trajectory {eq}`sdstraj` is - -```{math} -:label: sdslinmodpath - -x_0, \quad -a x_0 + b, \quad -a^2 x_0 + a b + b, \quad \text{etc.} -``` - -Continuing in this way, and using our knowledge of {doc}`geometric series `, we find that, for any $t \geq 0$, - -```{math} -:label: sdslinmod - -x_t = a^t x_0 + b \frac{1 - a^t}{1 - a} -``` - -This is about all we need to know about the linear model. - -We have an exact expression for $x_t$ for all $t$ and hence a full -understanding of the dynamics. - -Notice in particular that $|a| < 1$, then, by {eq}`sdslinmod`, we have - -```{math} -:label: sdslinmodc - -x_t \to \frac{b}{1 - a} \text{ as } t \to \infty -``` - -regardless of $x_0$ - -This is an example of what is called global stability, a topic we return to -below. - -### Example: A Nonlinear Model - -In the linear example above, we obtained an exact analytical expression for $x_t$ -in terms of arbitrary $t$ and $x_0$. - -This made analysis of dynamics very easy. - -When models are nonlinear, however, the situation can be quite different. - -For example, recall how we [previously studied](https://python-programming.quantecon.org/python_oop.html#example-the-solow-growth-model) the law of motion for the Solow growth model, a simplified version of which is - -```{math} -:label: solow_lom2 - -k_{t+1} = s z k_t^{\alpha} + (1 - \delta) k_t -``` - -Here $k$ is capital stock and $s, z, \alpha, \delta$ are positive -parameters with $0 < \alpha, \delta < 1$. - -If you try to iterate like we did in {eq}`sdslinmodpath`, you will find that -the algebra gets messy quickly. - -Analyzing the dynamics of this model requires a different method (see below). - -### Stability - -A **steady state** of the difference equation $x_{t+1} = g(x_t)$ is a -point $x^*$ in $S$ such that $x^* = g(x^*)$. - -In other words, $x^*$ is a **fixed point** of the function $g$ in -$S$. - -For example, for the linear model $x_{t+1} = a x_t + b$, you can use the -definition to check that - -* $x^* := b/(1-a)$ is a steady state whenever $a \not= 1$. -* if $a = 1$ and $b=0$, then every $x \in \mathbb R$ is a - steady state. -* if $a = 1$ and $b \not= 0$, then the linear model has no steady - state in $\mathbb R$. - -A steady state $x^*$ of $x_{t+1} = g(x_t)$ is called -**globally stable** if, for all $x_0 \in S$, - -$$ -x_t = g^t(x_0) \to x^* \text{ as } t \to \infty -$$ - -For example, in the linear model $x_{t+1} = a x_t + b$ with $a -\not= 1$, the steady state $x^*$ - -* is globally stable if $|a| < 1$ and -* fails to be globally stable otherwise. - -This follows directly from {eq}`sdslinmod`. - -A steady state $x^*$ of $x_{t+1} = g(x_t)$ is called -**locally stable** if there exists an $\epsilon > 0$ such that - -$$ -| x_0 - x^* | < \epsilon -\; \implies \; -x_t = g^t(x_0) \to x^* \text{ as } t \to \infty -$$ - -Obviously every globally stable steady state is also locally stable. - -We will see examples below where the converse is not true. - -## Graphical Analysis - -As we saw above, analyzing the dynamics for nonlinear models is nontrivial. - -There is no single way to tackle all nonlinear models. - -However, there is one technique for one-dimensional models that provides a -great deal of intuition. - -This is a graphical approach based on **45 degree diagrams**. - -Let's look at an example: the Solow model with dynamics given in {eq}`solow_lom2`. - -We begin with some plotting code that you can ignore at first reading. - -The function of the code is to produce 45 degree diagrams and time series -plots. - -```{code-cell} ipython ---- -tags: [output_scroll] ---- -def subplots(fs): - "Custom subplots with axes throught the origin" - fig, ax = plt.subplots(figsize=fs) - - # Set the axes through the origin - for spine in ['left', 'bottom']: - ax.spines[spine].set_position('zero') - ax.spines[spine].set_color('green') - for spine in ['right', 'top']: - ax.spines[spine].set_color('none') - - return fig, ax - - -def plot45(g, xmin, xmax, x0, num_arrows=6, var='x'): - - xgrid = np.linspace(xmin, xmax, 200) - - fig, ax = subplots((6.5, 6)) - ax.set_xlim(xmin, xmax) - ax.set_ylim(xmin, xmax) - - hw = (xmax - xmin) * 0.01 - hl = 2 * hw - arrow_args = dict(fc="k", ec="k", head_width=hw, - length_includes_head=True, lw=1, - alpha=0.6, head_length=hl) - - ax.plot(xgrid, g(xgrid), 'b-', lw=2, alpha=0.6, label='g') - ax.plot(xgrid, xgrid, 'k-', lw=1, alpha=0.7, label='45') - - x = x0 - xticks = [xmin] - xtick_labels = [xmin] - - for i in range(num_arrows): - if i == 0: - ax.arrow(x, 0.0, 0.0, g(x), **arrow_args) # x, y, dx, dy - else: - ax.arrow(x, x, 0.0, g(x) - x, **arrow_args) - ax.plot((x, x), (0, x), 'k', ls='dotted') - - ax.arrow(x, g(x), g(x) - x, 0, **arrow_args) - xticks.append(x) - xtick_labels.append(r'${}_{}$'.format(var, str(i))) - - x = g(x) - xticks.append(x) - xtick_labels.append(r'${}_{}$'.format(var, str(i+1))) - ax.plot((x, x), (0, x), 'k', ls='dotted') - - xticks.append(xmax) - xtick_labels.append(xmax) - ax.set_xticks(xticks) - ax.set_yticks(xticks) - ax.set_xticklabels(xtick_labels) - ax.set_yticklabels(xtick_labels) - - bbox = (0., 1.04, 1., .104) - legend_args = {'bbox_to_anchor': bbox, 'loc': 'upper right'} - - ax.legend(ncol=2, frameon=False, **legend_args, fontsize=14) - plt.show() - -def ts_plot(g, xmin, xmax, x0, ts_length=6, var='x'): - fig, ax = subplots((7, 5.5)) - ax.set_ylim(xmin, xmax) - ax.set_xlabel(r'$t$', fontsize=14) - ax.set_ylabel(r'${}_t$'.format(var), fontsize=14) - x = np.empty(ts_length) - x[0] = x0 - for t in range(ts_length-1): - x[t+1] = g(x[t]) - ax.plot(range(ts_length), - x, - 'bo-', - alpha=0.6, - lw=2, - label=r'${}_t$'.format(var)) - ax.legend(loc='best', fontsize=14) - ax.set_xticks(range(ts_length)) - plt.show() -``` - -Let's create a 45 degree diagram for the Solow model with a fixed set of -parameters - -```{code-cell} ipython -A, s, alpha, delta = 2, 0.3, 0.3, 0.4 -``` - -Here's the update function corresponding to the model. - -```{code-cell} ipython -def g(k): - return A * s * k**alpha + (1 - delta) * k -``` - -Here is the 45 degree plot. - -```{code-cell} ipython -xmin, xmax = 0, 4 # Suitable plotting region. - -plot45(g, xmin, xmax, 0, num_arrows=0) -``` - -The plot shows the function $g$ and the 45 degree line. - -Think of $k_t$ as a value on the horizontal axis. - -To calculate $k_{t+1}$, we can use the graph of $g$ to see its -value on the vertical axis. - -Clearly, - -* If $g$ lies above the 45 degree line at this point, then we have $k_{t+1} > k_t$. -* If $g$ lies below the 45 degree line at this point, then we have $k_{t+1} < k_t$. -* If $g$ hits the 45 degree line at this point, then we have $k_{t+1} = k_t$, so $k_t$ is a steady state. - -For the Solow model, there are two steady states when $S = \mathbb R_+ = -[0, \infty)$. - -* the origin $k=0$ -* the unique positive number such that $k = s z k^{\alpha} + (1 - \delta) k$. - -By using some algebra, we can show that in the second case, the steady state is - -$$ -k^* = \left( \frac{sz}{\delta} \right)^{1/(1-\alpha)} -$$ - -### Trajectories - -By the preceding discussion, in regions where $g$ lies above the 45 degree line, we know that the trajectory is increasing. - -The next figure traces out a trajectory in such a region so we can see this more clearly. - -The initial condition is $k_0 = 0.25$. - -```{code-cell} ipython -k0 = 0.25 - -plot45(g, xmin, xmax, k0, num_arrows=5, var='k') -``` - -We can plot the time series of capital corresponding to the figure above as -follows: - -```{code-cell} ipython -ts_plot(g, xmin, xmax, k0, var='k') -``` - -Here's a somewhat longer view: - -```{code-cell} ipython -ts_plot(g, xmin, xmax, k0, ts_length=20, var='k') -``` - -When capital stock is higher than the unique positive steady state, we see that -it declines: - -```{code-cell} ipython -k0 = 2.95 - -plot45(g, xmin, xmax, k0, num_arrows=5, var='k') -``` - -Here is the time series: - -```{code-cell} ipython -ts_plot(g, xmin, xmax, k0, var='k') -``` - -### Complex Dynamics - -The Solow model is nonlinear but still generates very regular dynamics. - -One model that generates irregular dynamics is the **quadratic map** - -$$ -g(x) = 4 x (1 - x), -\qquad x \in [0, 1] -$$ - -Let's have a look at the 45 degree diagram. - -```{code-cell} ipython -xmin, xmax = 0, 1 -g = lambda x: 4 * x * (1 - x) - -x0 = 0.3 -plot45(g, xmin, xmax, x0, num_arrows=0) -``` - -Now let's look at a typical trajectory. - -```{code-cell} ipython -plot45(g, xmin, xmax, x0, num_arrows=6) -``` - -Notice how irregular it is. - -Here is the corresponding time series plot. - -```{code-cell} ipython -ts_plot(g, xmin, xmax, x0, ts_length=6) -``` - -The irregularity is even clearer over a longer time horizon: - -```{code-cell} ipython -ts_plot(g, xmin, xmax, x0, ts_length=20) -``` - -## Exercises - -```{exercise} -:label: sd_ex1 - -Consider again the linear model $x_{t+1} = a x_t + b$ with $a -\not=1$. - -The unique steady state is $b / (1 - a)$. - -The steady state is globally stable if $|a| < 1$. - -Try to illustrate this graphically by looking at a range of initial conditions. - -What differences do you notice in the cases $a \in (-1, 0)$ and $a -\in (0, 1)$? - -Use $a=0.5$ and then $a=-0.5$ and study the trajectories - -Set $b=1$ throughout. -``` - -```{solution-start} sd_ex1 -:class: dropdown -``` - -We will start with the case $a=0.5$. - -Let's set up the model and plotting region: - -```{code-cell} ipython -a, b = 0.5, 1 -xmin, xmax = -1, 3 -g = lambda x: a * x + b -``` - -Now let's plot a trajectory: - -```{code-cell} ipython -x0 = -0.5 -plot45(g, xmin, xmax, x0, num_arrows=5) -``` - -Here is the corresponding time series, which converges towards the steady -state. - -```{code-cell} ipython -ts_plot(g, xmin, xmax, x0, ts_length=10) -``` - -Now let's try $a=-0.5$ and see what differences we observe. - -Let's set up the model and plotting region: - -```{code-cell} ipython -a, b = -0.5, 1 -xmin, xmax = -1, 3 -g = lambda x: a * x + b -``` - -Now let's plot a trajectory: - -```{code-cell} ipython -x0 = -0.5 -plot45(g, xmin, xmax, x0, num_arrows=5) -``` - -Here is the corresponding time series, which converges towards the steady -state. - -```{code-cell} ipython -ts_plot(g, xmin, xmax, x0, ts_length=10) -``` - -Once again, we have convergence to the steady state but the nature of -convergence differs. - -In particular, the time series jumps from above the steady state to below it -and back again. - -In the current context, the series is said to exhibit **damped oscillations**. - -```{solution-end} -```