{
"cells": [
{
"cell_type": "markdown",
"id": "cd1dc6df",
"metadata": {},
"source": [
"# Applied Linear Algebra\n",
"\n",
"**Prerequisites**\n",
"\n",
"- [Introduction to Numpy](https://datascience.quantecon.org/numpy_arrays.html) \n",
"\n",
"\n",
"**Outcomes**\n",
"\n",
"- Refresh some important linear algebra concepts \n",
"- Apply concepts to understanding unemployment and pricing portfolios \n",
"- Use `numpy` to do linear algebra operations "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a5f3dd1e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Uncomment following line to install on colab\n",
"#! pip install "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "885970fa",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# import numpy to prepare for code below\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"id": "35288dfd",
"metadata": {},
"source": [
"## Vectors and Matrices"
]
},
{
"cell_type": "markdown",
"id": "569bf5d3",
"metadata": {},
"source": [
"### Vectors\n",
"\n",
"A (N-element) vector is $ N $ numbers stored together.\n",
"\n",
"We typically write a vector as $ x = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\dots \\\\ x_N \\end{bmatrix} $.\n",
"\n",
"In numpy terms, a vector is a 1-dimensional array.\n",
"\n",
"We often think of 2-element vectors as directional lines in the XY axes.\n",
"\n",
"This image, from the [QuantEcon Python lecture](https://python.quantecon.org/linear_algebra.html)\n",
"is an example of what this might look like for the vectors `(-4, 3.5)`, `(-3, 3)`, and `(2, 4)`.\n",
"\n",
"![https://datascience.quantecon.org/_static/vector.png](https://datascience.quantecon.org/_static/vector.png)\n",
"\n",
" \n",
"In a previous lecture, we saw some types of operations that can be done on\n",
"vectors, such as"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e8e200bf",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x = np.array([1, 2, 3])\n",
"y = np.array([4, 5, 6])"
]
},
{
"cell_type": "markdown",
"id": "db15b08c",
"metadata": {},
"source": [
"**Element-wise operations**: Let $ z = x ? y $ for some operation $ ? $, one of\n",
"the standard *binary* operations ($ +, -, \\times, \\div $). Then we can write\n",
"$ z = \\begin{bmatrix} x_1 ? y_1 & x_2 ? y_2 \\end{bmatrix} $. Element-wise operations require\n",
"that $ x $ and $ y $ have the same size."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f789e061",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Element-wise Addition\", x + y)\n",
"print(\"Element-wise Subtraction\", x - y)\n",
"print(\"Element-wise Multiplication\", x * y)\n",
"print(\"Element-wise Division\", x / y)"
]
},
{
"cell_type": "markdown",
"id": "c2f81864",
"metadata": {},
"source": [
"**Scalar operations**: Let $ w = a ? x $ for some operation $ ? $, one of the\n",
"standard *binary* operations ($ +, -, \\times, \\div $). Then we can write\n",
"$ w = \\begin{bmatrix} a ? x_1 & a ? x_2 \\end{bmatrix} $."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a22c391c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Scalar Addition\", 3 + x)\n",
"print(\"Scalar Subtraction\", 3 - x)\n",
"print(\"Scalar Multiplication\", 3 * x)\n",
"print(\"Scalar Division\", 3 / x)"
]
},
{
"cell_type": "markdown",
"id": "2ffed637",
"metadata": {},
"source": [
"Another operation very frequently used in data science is the **dot product**.\n",
"\n",
"The dot between $ x $ and $ y $ is written $ x \\cdot y $ and is\n",
"equal to $ \\sum_{i=1}^N x_i y_i $."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "88b12bea",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Dot product\", np.dot(x, y))"
]
},
{
"cell_type": "markdown",
"id": "f299548e",
"metadata": {},
"source": [
"We can also use `@` to denote dot products (and matrix multiplication which we’ll see soon!)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19536929",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Dot product with @\", x @ y)"
]
},
{
"cell_type": "markdown",
"id": "b08fffc4",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"See exercise 1 in the [exercise list](#ex3-3)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7461607a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"nA = 100\n",
"nB = 50\n",
"nassets = np.array([nA, nB])\n",
"\n",
"i = 0.05\n",
"durationA = 6\n",
"durationB = 4\n",
"\n",
"# Do your computations here\n",
"\n",
"# Compute price\n",
"\n",
"# uncomment below to see a message!\n",
"# if condition:\n",
"# print(\"Alice can retire\")\n",
"# else:\n",
"# print(\"Alice cannot retire yet\")"
]
},
{
"cell_type": "markdown",
"id": "0593f125",
"metadata": {},
"source": [
"### Matrices\n",
"\n",
"An $ N \\times M $ matrix can be thought of as a collection of M\n",
"N-element vectors stacked side-by-side as columns.\n",
"\n",
"We write a matrix as\n",
"\n",
"$$\n",
"\\begin{bmatrix} x_{11} & x_{12} & \\dots & x_{1M} \\\\\n",
" x_{21} & \\dots & \\dots & x_{2M} \\\\\n",
" \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
" x_{N1} & x_{N2} & \\dots & x_{NM}\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"In numpy terms, a matrix is a 2-dimensional array.\n",
"\n",
"We can create a matrix by passing a list of lists to the `np.array` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f79b1e1f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x = np.array([[1, 2, 3], [4, 5, 6]])\n",
"y = np.ones((2, 3))\n",
"z = np.array([[1, 2], [3, 4], [5, 6]])"
]
},
{
"cell_type": "markdown",
"id": "a493bf19",
"metadata": {},
"source": [
"We can perform element-wise and scalar operations as we did with vectors. In fact, we can do\n",
"these two operations on arrays of any dimension."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5f2b020",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Element-wise Addition\\n\", x + y)\n",
"print(\"Element-wise Subtraction\\n\", x - y)\n",
"print(\"Element-wise Multiplication\\n\", x * y)\n",
"print(\"Element-wise Division\\n\", x / y)\n",
"\n",
"print(\"Scalar Addition\\n\", 3 + x)\n",
"print(\"Scalar Subtraction\\n\", 3 - x)\n",
"print(\"Scalar Multiplication\\n\", 3 * x)\n",
"print(\"Scalar Division\\n\", 3 / x)"
]
},
{
"cell_type": "markdown",
"id": "cf481dd5",
"metadata": {},
"source": [
"Similar to how we combine vectors with a dot product, matrices can do what we’ll call *matrix\n",
"multiplication*.\n",
"\n",
"Matrix multiplication is effectively a generalization of dot products.\n",
"\n",
"**Matrix multiplication**: Let $ v = x \\cdot y $ then we can write\n",
"$ v_{ij} = \\sum_{k=1}^N x_{ik} y_{kj} $ where $ x_{ij} $ is notation that denotes the\n",
"element found in the ith row and jth column of the matrix $ x $.\n",
"\n",
"The image below from [Wikipedia](https://commons.wikimedia.org/wiki/File:Matrix_multiplication_diagram.svg),\n",
"by Bilou, shows how matrix multiplication simplifies to a series of dot products:\n",
"\n",
"![https://datascience.quantecon.org/_static/mat_mult_wiki_bilou.png](https://datascience.quantecon.org/_static/mat_mult_wiki_bilou.png)\n",
"\n",
" \n",
"After looking at the math and image above, you might have realized that matrix\n",
"multiplication requires very specific matrix shapes!\n",
"\n",
"For two matrices $ x, y $ to be multiplied, $ x $\n",
"must have the same number of columns as $ y $ has rows.\n",
"\n",
"Formally, we require that for some integer numbers, $ M, N, $ and $ K $\n",
"that if $ x $ is $ N \\times M $ then $ y $ must be $ M \\times\n",
"K $.\n",
"\n",
"If we think of a vector as a $ 1 \\times M $ or $ M \\times 1 $ matrix, we can even do\n",
"matrix multiplication between a matrix and a vector!\n",
"\n",
"Let’s see some examples of this."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13b0db56",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x1 = np.reshape(np.arange(6), (3, 2))\n",
"x2 = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])\n",
"x3 = np.array([[2, 5, 2], [1, 2, 1]])\n",
"x4 = np.ones((2, 3))\n",
"\n",
"y1 = np.array([1, 2, 3])\n",
"y2 = np.array([0.5, 0.5])"
]
},
{
"cell_type": "markdown",
"id": "149c8d23",
"metadata": {},
"source": [
"Numpy allows us to do matrix multiplication in three ways."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a18b1ca0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Using the matmul function for two matrices\")\n",
"print(np.matmul(x1, x4))\n",
"print(\"Using the dot function for two matrices\")\n",
"print(np.dot(x1, x4))\n",
"print(\"Using @ for two matrices\")\n",
"print(x1 @ x4)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f177ea09",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(\"Using the matmul function for vec and mat\")\n",
"print(np.matmul(y1, x1))\n",
"print(\"Using the dot function for vec and mat\")\n",
"print(np.dot(y1, x1))\n",
"print(\"Using @ for vec and mat\")\n",
"print(y1 @ x1)"
]
},
{
"cell_type": "markdown",
"id": "af731a54",
"metadata": {},
"source": [
"Despite our options, we stick to using `@` because\n",
"it is simplest to read and write."
]
},
{
"cell_type": "markdown",
"id": "09c51bfe",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"See exercise 2 in the [exercise list](#ex3-3)."
]
},
{
"cell_type": "markdown",
"id": "6f690331",
"metadata": {},
"source": [
"### Other Linear Algebra Concepts"
]
},
{
"cell_type": "markdown",
"id": "f85a946b",
"metadata": {},
"source": [
"#### Transpose\n",
"\n",
"A matrix transpose is an operation that flips all elements of a matrix along the diagonal.\n",
"\n",
"More formally, the $ (i, j) $ element of $ x $ becomes the $ (j, i) $ element of\n",
"$ x^T $.\n",
"\n",
"In particular, let $ x $ be given by\n",
"\n",
"$$\n",
"x = \\begin{bmatrix} 1 & 2 & 3 \\\\\n",
" 4 & 5 & 6 \\\\\n",
" 7 & 8 & 9 \\\\\n",
" \\end{bmatrix}\n",
"$$\n",
"\n",
"then $ x $ transpose, written as $ x' $, is given by\n",
"\n",
"$$\n",
"x = \\begin{bmatrix} 1 & 4 & 7 \\\\\n",
" 2 & 5 & 8 \\\\\n",
" 3 & 6 & 9 \\\\\n",
" \\end{bmatrix}\n",
"$$\n",
"\n",
"In Python, we do this by"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1d2d390d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
"\n",
"print(\"x transpose is\")\n",
"print(x.transpose())"
]
},
{
"cell_type": "markdown",
"id": "c1459f36",
"metadata": {},
"source": [
"#### Identity Matrix\n",
"\n",
"In linear algebra, one particular matrix acts very similarly to how 1 behaves for scalar numbers.\n",
"\n",
"This matrix is known as the *identity matrix* and is given by\n",
"\n",
"$$\n",
"I = \\begin{bmatrix} 1 & 0 & 0 & \\dots & 0 \\\\\n",
" 0 & 1 & 0 & \\dots & 0 \\\\\n",
" \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
" 0 & 0 & 0 & \\dots & 1\n",
" \\end{bmatrix}\n",
"$$\n",
"\n",
"As seen above, it has 1s on the diagonal and 0s everywhere else.\n",
"\n",
"When we multiply any matrix or vector by the identity matrix, we get the original matrix or vector\n",
"back!\n",
"\n",
"Let’s see some examples."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5d3b349c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"I = np.eye(3)\n",
"x = np.reshape(np.arange(9), (3, 3))\n",
"y = np.array([1, 2, 3])\n",
"\n",
"print(\"I @ x\", \"\\n\", I @ x)\n",
"print(\"x @ I\", \"\\n\", x @ I)\n",
"print(\"I @ y\", \"\\n\", I @ y)\n",
"print(\"y @ I\", \"\\n\", y @ I)"
]
},
{
"cell_type": "markdown",
"id": "1f14aba6",
"metadata": {},
"source": [
"#### Inverse\n",
"\n",
"If you recall, you learned in your primary education about solving equations for certain variables.\n",
"\n",
"For example, you might have been given the equation\n",
"\n",
"$$\n",
"3x + 7 = 16\n",
"$$\n",
"\n",
"and then asked to solve for $ x $.\n",
"\n",
"You probably did this by subtracting 7 and then dividing by 3.\n",
"\n",
"Now let’s write an equation that contains matrices and vectors.\n",
"\n",
"$$\n",
"\\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix} \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} = \\begin{bmatrix} 3 \\\\ 4 \\end{bmatrix}\n",
"$$\n",
"\n",
"How would we solve for $ x = \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} $?\n",
"\n",
"Unfortunately, there is no “matrix divide” operation that does the opposite of matrix multiplication.\n",
"\n",
"Instead, we first have to do what’s known as finding the inverse. We must multiply both sides by this inverse to solve.\n",
"\n",
"Consider some matrix $ A $.\n",
"\n",
"The inverse of $ A $, given by $ A^{-1} $, is a matrix such that $ A A^{-1} = I $\n",
"where $ I $ is our identity matrix.\n",
"\n",
"Notice in our equation above, if we can find the inverse of\n",
"$ \\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix} $ then we can multiply both sides by the inverse\n",
"to get\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix}^{-1}\\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix} \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} &= \\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix}^{-1}\\begin{bmatrix} 3 \\\\ 4 \\end{bmatrix} \\\\\n",
"I \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} &= \\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix}^{-1} \\begin{bmatrix} 3 \\\\ 4 \\end{bmatrix} \\\\\n",
" \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} &= \\begin{bmatrix} 1 & 2 \\\\ 3 & 1 \\end{bmatrix}^{-1} \\begin{bmatrix} 3 \\\\ 4 \\end{bmatrix}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Computing the inverse requires that a matrix be square and satisfy some other conditions\n",
"(non-singularity) that are beyond the scope of this lecture.\n",
"\n",
"We also skip the exact details of how this inverse is computed, but, if you are interested,\n",
"you can visit the\n",
"[QuantEcon Linear Algebra lecture](https://python.quantecon.org/linear_algebra.html)\n",
"for more details.\n",
"\n",
"We demonstrate how to compute the inverse with numpy below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2b6615a0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# This is a square (N x N) non-singular matrix\n",
"A = np.array([[1, 2, 0], [3, 1, 0], [0, 1, 2]])\n",
"\n",
"print(\"This is A inverse\")\n",
"\n",
"print(np.linalg.inv(A))\n",
"\n",
"print(\"Check that A @ A inverse is I\")\n",
"print(np.linalg.inv(A) @ A)"
]
},
{
"cell_type": "markdown",
"id": "cfa851a2",
"metadata": {},
"source": [
"## Portfolios\n",
"\n",
"In [control flow](https://datascience.quantecon.org/../python_fundamentals/control_flow.html), we learned to value a stream of payoffs from a single\n",
"asset.\n",
"\n",
"In this section, we generalize this to value a portfolio of multiple assets, or an asset\n",
"that has easily separable components.\n",
"\n",
"Vectors and inner products give us a convenient way to organize and calculate these payoffs."
]
},
{
"cell_type": "markdown",
"id": "a1081002",
"metadata": {},
"source": [
"### Static Payoffs\n",
"\n",
"As an example, consider a portfolio with 4 units of asset A, 2.5 units of asset B, and 8 units of\n",
"asset C.\n",
"\n",
"At a particular point in time, the assets pay $ 3 $/unit of asset A, $ 5 $/unit of B, and\n",
"$ 1.10 $/unit of C.\n",
"\n",
"First, calculate the value of this portfolio directly with a sum."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3347060e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"4.0 * 3.0 + 2.5 * 5.0 + 8 * 1.1"
]
},
{
"cell_type": "markdown",
"id": "cb59a376",
"metadata": {},
"source": [
"We can make this more convenient and general by using arrays for accounting, and then sum then in a\n",
"loop."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "03983615",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"x = np.array([4.0, 2.5, 8.0]) # portfolio units\n",
"y = np.array([3.0, 5.0, 1.1]) # payoffs\n",
"n = len(x)\n",
"p = 0.0\n",
"for i in range(n): # i.e. 0, 1, 2\n",
" p = p + x[i] * y[i]\n",
"\n",
"p"
]
},
{
"cell_type": "markdown",
"id": "c0f25479",
"metadata": {},
"source": [
"The above would have worked with `x` and `y` as `list` rather than `np.array`.\n",
"\n",
"Note that the general pattern above is the sum.\n",
"\n",
"$$\n",
"p = \\sum_{i=0}^{n-1} x_i y_i = x \\cdot y\n",
"$$\n",
"\n",
"This is an inner product as implemented by the `np.dot` function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e19ce2be",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"np.dot(x, y)"
]
},
{
"cell_type": "markdown",
"id": "e024506f",
"metadata": {},
"source": [
"This approach allows us to simultaneously price different portfolios by stacking them in a matrix and using the dot product."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0242c611",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"y = np.array([3.0, 5.0, 1.1]) # payoffs\n",
"x1 = np.array([4.0, 2.5, 8.0]) # portfolio 1\n",
"x2 = np.array([2.0, 1.5, 0.0]) # portfolio 2\n",
"X = np.array((x1, x2))\n",
"\n",
"# calculate with inner products\n",
"p1 = np.dot(X[0,:], y)\n",
"p2 = np.dot(X[1,:], y)\n",
"print(\"Calculating separately\")\n",
"print([p1, p2])\n",
"\n",
"# or with a matrix multiplication\n",
"print(\"Calculating with matrices\")\n",
"P = X @ y\n",
"print(P)"
]
},
{
"cell_type": "markdown",
"id": "a594e59a",
"metadata": {},
"source": [
"### NPV of a Portfolio\n",
"\n",
"If a set of assets has payoffs over time, we can calculate the NPV of that portfolio in a similar way to the calculation in\n",
"[npv](https://datascience.quantecon.org/../python_fundamentals/control_flow.html#npv).\n",
"\n",
"First, consider an example with an asset with claims to multiple streams of payoffs which are easily\n",
"separated.\n",
"\n",
"You are considering purchasing an oilfield with 2 oil wells, named `A` and `B` where\n",
"\n",
"- Both oilfields have a finite lifetime of 20 years. \n",
"- In oilfield `A`, you can extract 5 units in the first year, and production in each subsequent year\n",
" decreases by $ 20\\% $ of the previous year so that\n",
" $ x^A_0 = 5, x^A_1 = 0.8 \\times 5, x^A_2 = 0.8^2 \\times 5, \\ldots $ \n",
"- In oilfield `B`, you can extract 2 units in the first year, but production only drops by\n",
" $ 10\\% $ each year (i.e. $ x^B_0 = 2, x^B_1 = 0.9 \\times 2, x^B_2 = 0.9^2 \\times 2, \\ldots $ \n",
"- Future cash flows are discounted at a rate of $ r = 0.05 $ each year. \n",
"- The price for oil in both wells are normalized as $ p_A = p_B = 1 $. \n",
"\n",
"\n",
"These traits can be separated so that the price you would be willing to pay is the sum of the two, where\n",
"we define $ \\gamma_A = 0.8, \\gamma_B = 0.9 $.\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"V_A &= \\sum_{t=0}^{T-1} \\left(\\frac{1}{1 + r}\\right)^t p_A y^A_t = \\sum_{t=0}^{T-1} \\left(\\frac{1}{1 + r}\\right)^t (p_A \\, x_{A0}\\, \\gamma_A^t)\\\\\n",
"V_B &= \\sum_{t=0}^{T-1} \\left(\\frac{1}{1 + r}\\right)^t p_B y^B_t = \\sum_{t=0}^{T-1} \\left(\\frac{1}{1 + r}\\right)^t (p_B \\, x_{B0}\\, \\gamma_B^t)\\\\\n",
"V &= V_A + V_B\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Let’s compute the value of each of these assets using the dot product.\n",
"\n",
"The first question to ask yourself is: “For which two vectors should I compute the dot product?”\n",
"\n",
"It turns out that this depends on which two vectors you’d like to create.\n",
"\n",
"One reasonable choice is presented in the code below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f46b97f4",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Depreciation of production rates\n",
"gamma_A = 0.80\n",
"gamma_B = 0.90\n",
"\n",
"# Interest rate discounting\n",
"r = 0.05\n",
"discount = np.array([(1 / (1+r))**t for t in range(20)])\n",
"\n",
"# Let's first create arrays that have the production of each oilfield\n",
"oil_A = 5 * np.array([gamma_A**t for t in range(20)])\n",
"oil_B = 2 * np.array([gamma_B**t for t in range(20)])\n",
"oilfields = np.array([oil_A, oil_B])\n",
"\n",
"# Use matrix multiplication to get discounted sum of oilfield values and then sum\n",
"# the two values\n",
"Vs = oilfields @ discount\n",
"\n",
"print(f\"The npv of oilfields is {Vs.sum()}\")"
]
},
{
"cell_type": "markdown",
"id": "9b473c82",
"metadata": {},
"source": [
"Now consider the approximation where instead of the oilfields having a finite lifetime of 20 years,\n",
"we let them produce forever, i.e. $ T = \\infty $.\n",
"\n",
"With a little algebra,\n",
"\n",
"$$\n",
"V_A = p_A \\sum_{t=0}^{\\infty}\\left(\\frac{1}{1 + r}\\right)^t (x_{A0} \\gamma_A^t) = x_{A0}\\sum_{t=0}^{\\infty}\\left(\\frac{\\gamma_A}{1 + r}\\right)^t\n",
"$$\n",
"\n",
"And, using the infinite sum formula from [Control Flow](https://datascience.quantecon.org/../python_fundamentals/control_flow.html) (i.e. $ \\sum_{t=0}^{\\infty}\\beta^t = (1 - \\beta)^{-1} $)\n",
"\n",
"$$\n",
"= \\frac{p_A x_{A0}}{1 - \\left(\\gamma_A\\frac{1}{1 + r} \\right)}\n",
"$$\n",
"\n",
"The $ V_B $ is defined symmetrically.\n",
"\n",
"How different is this infinite horizon approximation from the $ T = 20 $ version, and why?\n",
"\n",
"Now, let’s compute the $ T = \\infty $ version of the net present value and make a graph to help\n",
"us see how many periods are needed to approach the infinite horizon value."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "77869ba3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Depreciation of production rates\n",
"gamma_A = 0.80\n",
"gamma_B = 0.90\n",
"\n",
"# Interest rate discounting\n",
"r = 0.05\n",
"\n",
"\n",
"def infhor_NPV_oilfield(starting_output, gamma, r):\n",
" beta = gamma / (1 + r)\n",
" return starting_output / (1 - beta)\n",
"\n",
"\n",
"def compute_NPV_oilfield(starting_output, gamma, r, T):\n",
" outputs = starting_output * np.array([gamma**t for t in range(T)])\n",
" discount = np.array([(1 / (1+r))**t for t in range(T)])\n",
"\n",
" npv = np.dot(outputs, discount)\n",
"\n",
" return npv\n",
"\n",
"Ts = np.arange(2, 75)\n",
"\n",
"NPVs_A = np.array([compute_NPV_oilfield(5, gamma_A, r, t) for t in Ts])\n",
"NPVs_B = np.array([compute_NPV_oilfield(2, gamma_B, r, t) for t in Ts])\n",
"\n",
"NPVs_T = NPVs_A + NPVs_B\n",
"NPV_oo = infhor_NPV_oilfield(5, gamma_A, r) + infhor_NPV_oilfield(2, gamma_B, r)\n",
"\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.set_title(\"NPV with Varying T\")\n",
"ax.set_ylabel(\"NPV\")\n",
"\n",
"ax.plot(Ts, NPVs_A + NPVs_B)\n",
"ax.hlines(NPV_oo, Ts[0], Ts[-1], color=\"k\", linestyle=\"--\") # Plot infinite horizon value\n",
"\n",
"ax.spines[\"right\"].set_visible(False)\n",
"ax.spines[\"top\"].set_visible(False)"
]
},
{
"cell_type": "markdown",
"id": "8003080b",
"metadata": {},
"source": [
"It is also worth noting that the computation of the infinite horizon net present value can be\n",
"simplified even further by using matrix multiplication. That is, the formula given above is\n",
"equivalent to\n",
"\n",
"$$\n",
"V = \\begin{bmatrix}p_A & p_B \\end{bmatrix} \\cdot \\sum_{t=0}^{\\infty} \\left(\\left(\\frac{1}{1 + r}\\right)^t \\begin{bmatrix} \\gamma_A & 0 \\\\ 0 & \\gamma_B \\end{bmatrix}^t \\cdot x_0\\right)\n",
"$$\n",
"\n",
"and where $ x_0 = \\begin{bmatrix} x_{A0} \\\\ x_{B0} \\end{bmatrix} $.\n",
"\n",
"We recognize that this equation is of the form\n",
"\n",
"$$\n",
"V = G \\sum_{t=0}^{\\infty} \\left(\\frac{1}{1 + r}\\right)^t A^t x_0\n",
"$$\n",
"\n",
"Without proof, and given important assumptions on $ \\frac{1}{1 + r} $ and $ A $, this\n",
"equation reduces to\n",
"\n",
"\n",
"\n",
"$$\n",
"V = G \\left(I - \\frac{1}{1+r} A\\right)^{-1} x_0 \\tag{1}\n",
"$$\n",
"\n",
"Using the matrix inverse, where `I` is the identity matrix."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "451f84f5",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"p_A = 1.0\n",
"p_B = 1.0\n",
"G = np.array([p_A, p_B])\n",
"\n",
"r = 0.05\n",
"beta = 1 / (1 + r)\n",
"\n",
"gamma_A = 0.80\n",
"gamma_B = 0.90\n",
"A = np.array([[gamma_A, 0], [0, gamma_B]])\n",
"\n",
"x_0 = np.array([5, 2])\n",
"\n",
"# Compute with matrix formula\n",
"NPV_mf = G @ np.linalg.inv(np.eye(2) - beta*A) @ x_0\n",
"\n",
"print(NPV_mf)"
]
},
{
"cell_type": "markdown",
"id": "ffd20eac",
"metadata": {},
"source": [
"Note: While our matrix above was very simple, this approach works for much more\n",
"complicated `A` matrices as long as we can write $ x_t $ using $ A $ and $ x_0 $ as\n",
"$ x_t = A^t x_0 $ (For an advanced description of this topic, adding randomness, read about\n",
"linear state-space models with Python [https://python.quantecon.org/linear_models.html](https://python.quantecon.org/linear_models.html))."
]
},
{
"cell_type": "markdown",
"id": "2fd2fda2",
"metadata": {},
"source": [
"### Unemployment Dynamics\n",
"\n",
"Consider an economy where in any given year, $ \\alpha = 5\\% $ of workers lose their jobs and\n",
"$ \\phi = 10\\% $ of unemployed workers find jobs.\n",
"\n",
"Define the vector $ x_0 = \\begin{bmatrix} 900,000 & 100,000 \\end{bmatrix} $ as the number of\n",
"employed and unemployed workers (respectively) at time $ 0 $ in the economy.\n",
"\n",
"Our goal is to determine the dynamics of unemployment in this economy.\n",
"\n",
"First, let’s define the matrix.\n",
"\n",
"$$\n",
"A = \\begin{bmatrix} 1 - \\alpha & \\alpha \\\\ \\phi & 1 - \\phi \\end{bmatrix}\n",
"$$\n",
"\n",
"Note that with this definition, we can describe the evolution of employment and unemployment\n",
"from $ x_0 $ to $ x_1 $ using linear algebra.\n",
"\n",
"$$\n",
"x_1 = \\begin{bmatrix} (1 - \\alpha) 900,000 + \\phi 100,000 \\\\ \\alpha 900,000 + (1-\\phi) 100,000\\end{bmatrix} = A' x_0\n",
"$$\n",
"\n",
"However, since the transitions do not change over time, we can use this to describe the evolution\n",
"from any arbitrary time $ t $, so that\n",
"\n",
"$$\n",
"x_{t+1} = A' x_t\n",
"$$\n",
"\n",
"Let’s code up a python function that will let us track the evolution of unemployment over time."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6910801b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"phi = 0.1\n",
"alpha = 0.05\n",
"\n",
"x0 = np.array([900_000, 100_000])\n",
"\n",
"A = np.array([[1-alpha, alpha], [phi, 1-phi]])\n",
"\n",
"def simulate(x0, A, T=10):\n",
" \"\"\"\n",
" Simulate the dynamics of unemployment for T periods starting from x0\n",
" and using values of A for probabilities of moving between employment\n",
" and unemployment\n",
" \"\"\"\n",
" nX = x0.shape[0]\n",
" out = np.zeros((T, nX))\n",
" out[0, :] = x0\n",
"\n",
" for t in range(1, T):\n",
" out[t, :] = A.T @ out[t-1, :]\n",
"\n",
" return out"
]
},
{
"cell_type": "markdown",
"id": "1018267f",
"metadata": {},
"source": [
"Let’s use this function to plot unemployment and employment levels for 10 periods."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "823d5e0b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def plot_simulation(x0, A, T=100):\n",
" X = simulate(x0, A, T)\n",
" fig, ax = plt.subplots()\n",
" ax.plot(X[:, 0])\n",
" ax.plot(X[:, 1])\n",
" ax.set_xlabel(\"t\")\n",
" ax.legend([\"Employed\", \"Unemployed\"])\n",
" return ax\n",
"\n",
"plot_simulation(x0, A, 50)"
]
},
{
"cell_type": "markdown",
"id": "275ae9c4",
"metadata": {},
"source": [
"Notice that the levels of unemployed an employed workers seem to be heading to constant numbers.\n",
"\n",
"We refer to this phenomenon as *convergence* because the values appear to converge to a constant\n",
"number.\n",
"\n",
"Let’s check that the values are permanently converging."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4332964f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_simulation(x0, A, 5000)"
]
},
{
"cell_type": "markdown",
"id": "8f9d0e0f",
"metadata": {},
"source": [
"The convergence of this system is a property determined by the matrix $ A $.\n",
"\n",
"The long-run distribution of employed and unemployed workers is equal to the [eigenvector](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors)\n",
"of $ A' $ that corresponds with the eigenvalue equal to 1. (An eigenvector of $ A' $ is also known as a “left-eigenvector” of $ A $.)\n",
"\n",
"Let’s have numpy compute the eigenvalues and eigenvectors and compare the results to our simulated results above:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10802541",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"eigvals, eigvecs = np.linalg.eig(A.T)\n",
"for i in range(len(eigvals)):\n",
" if eigvals[i] == 1:\n",
" which_eig = i\n",
" break\n",
"\n",
"print(f\"We are looking for eigenvalue {which_eig}\")"
]
},
{
"cell_type": "markdown",
"id": "cac45c92",
"metadata": {},
"source": [
"Now let’s look at the corresponding eigenvector:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3a7af5d4",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"dist = eigvecs[:, which_eig]\n",
"\n",
"# need to divide by sum so it adds to 1\n",
"dist /= dist.sum()\n",
"\n",
"print(f\"The distribution of workers is given by {dist}\")"
]
},
{
"cell_type": "markdown",
"id": "9b5c9177",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"See exercise 3 in the [exercise list](#ex3-3).\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "064d6397",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"id": "bbf59667",
"metadata": {},
"source": [
"### Exercise 1\n",
"\n",
"Alice is a stock broker who owns two types of assets: A and B. She owns 100\n",
"units of asset A and 50 units of asset B. The current interest rate is 5%.\n",
"Each of the A assets have a remaining duration of 6 years and pay\n",
"\\$1500 each year, while each of the B assets have a remaining duration\n",
"of 4 years and pay \\$500 each year. Alice would like to retire if she\n",
"can sell her assets for more than \\$500,000. Use vector addition, scalar\n",
"multiplication, and dot products to determine whether she can retire.\n",
"\n",
"([back to text](#dir3-3-1))"
]
},
{
"cell_type": "markdown",
"id": "73f0637c",
"metadata": {},
"source": [
"### Exercise 2\n",
"\n",
"Which of the following operations will work and which will\n",
"create errors because of size issues?\n",
"\n",
"Test out your intuitions in the code cell below"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ab970686",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x1 @ x2\n",
"x2 @ x1\n",
"x2 @ x3\n",
"x3 @ x2\n",
"x1 @ x3\n",
"x4 @ y1\n",
"x4 @ y2\n",
"y1 @ x4\n",
"y2 @ x4"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12321cca",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# testing area"
]
},
{
"cell_type": "markdown",
"id": "36bb8c41",
"metadata": {},
"source": [
"([back to text](#dir3-3-2))"
]
},
{
"cell_type": "markdown",
"id": "4e26581b",
"metadata": {},
"source": [
"### Exercise 3\n",
"\n",
"Compare the distribution above to the final values of a long simulation.\n",
"\n",
"If you multiply the distribution by 1,000,000 (the number of workers), do you get (roughly) the same number as the simulation?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f9fe4a39",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# your code here"
]
},
{
"cell_type": "markdown",
"id": "add4ff53",
"metadata": {},
"source": [
"([back to text](#dir3-3-3))"
]
}
],
"metadata": {
"date": 1689807041.9201207,
"filename": "applied_linalg.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Applied Linear Algebra"
},
"nbformat": 4,
"nbformat_minor": 5
}