maths-numerical.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bf4e32e9",
   "metadata": {},
   "source": [
    "# Numerical Mathematics\n",
    "\n",
    "In this chapter, you'll learn about how to solve mathematical problems numerically.\n",
    "\n",
    "If you're running code from this chapter, remember you may need to install the packages. This chapter uses **numpy**, **scipy**, and **matplotlib** in particular."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e1f0ee",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "88cca34b",
   "metadata": {},
   "source": [
    "## Arrays\n",
    "\n",
    "**Numpy** arrays are the fundamental non-scalar object for numerical analysis in Python. Let's meet a vector of integers rendered as a **numpy** array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a61ecb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "a = np.array([0, 1, 2, 3], dtype=\"int64\")\n",
    "a"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1243113b",
   "metadata": {},
   "source": [
    "Arrays are very memory efficient and fast objects that you should use in preference to lists for any heavy duty numerical operation.\n",
    "\n",
    "To demonstrate this, let's do a time race between lists and arrays for squaring all elements of an array:\n",
    "\n",
    "Lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5246ee0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_list = range(1000)\n",
    "%timeit [i**2 for i in a_list]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8f8d01e0",
   "metadata": {},
   "source": [
    "Arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c40b6403",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1000)\n",
    "%timeit a**2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "26ce79d4",
   "metadata": {},
   "source": [
    "Using arrays was *two orders of magnitude** faster! Okay, so we should use arrays for numerical works. How do we make them? You can specify an array explicitly as we did above to create a vector. This manual approach works for other dimensions too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde5edcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])\n",
    "mat"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "96be90e2",
   "metadata": {},
   "source": [
    "To find out about matrix properties, we use `.shape`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00f62dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "mat.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d841b4bd",
   "metadata": {},
   "source": [
    "We already saw how `np.arange(start, stop, step)` produces a vector; `np.linspace(start, stop, number)` produces a vector of length `number` by equally dividing the space between `start` and `stop`.\n",
    "\n",
    "Three really useful arrays are `np.ones(shape)`, for example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ddb5e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones((3, 3))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "08bab38a",
   "metadata": {},
   "source": [
    "`np.diag()` for diagonal arrays,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "373af9c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.diag(np.array([1, 2, 3, 4]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a08e557a",
   "metadata": {},
   "source": [
    "and `np.zeros()` for empty arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f67a75a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros((2, 2))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0f01e3d7",
   "metadata": {},
   "source": [
    "Random numbers are supplied by `np.random.rand()` for a uniform distribution in [0, 1], and `np.random.randn()` for numbers drawn from a standard normal distribution.\n",
    "\n",
    "You can, of course, specify a function to create an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4f4e5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.fromfunction(lambda i, j: i**2 + j**2, (4, 5))\n",
    "c"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6914be94",
   "metadata": {},
   "source": [
    "To access values in an array, you can use all of the by-position slicing methods that you've seen already in data analysis and with lists. The figure gives an example of some common slicing operations:\n",
    "\n",
    "![Examples of array slices](https://scipy-lectures.org/_images/numpy_indexing.png)\n",
    "\n",
    "Arrays can also be sliced and diced based on boolean indexing, just like a dataframe.\n",
    "\n",
    "For example, using the array defined above, we can create a boolean array of true and false values from a condition such as `c > 6` and use that to only access some elements of an array (it doesn't have to be the same array, though it usually is):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c729b2b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "c[c > 6]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d3498beb",
   "metadata": {},
   "source": [
    "As with dataframes, arrays can be combined. The main command to remember is `np.concatenate`, which has an `axis` keyword option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f35351",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.eye(3)\n",
    "np.concatenate([x, x], axis=0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "280a5a19",
   "metadata": {},
   "source": [
    "Splitting is performed with `np.split(array, splits, axis=)`, for example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad1f449f",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.split(x, [3], axis=0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c84b1066",
   "metadata": {},
   "source": [
    "Aggregation operations are very similar to those found in dataframes: `x.sum(i)` to sum across the $i$th dimension of the array; similarly for standard deviation, and so on.\n",
    "\n",
    "As with dataframes, you can (and often should) specify the datatype of an array when you create it by passing a `dtype=` keyword, eg `c = np.array([1, 2, 3], dtype=float)`. To find out the data type of an array that already exists, use `c.dtype`.\n",
    "\n",
    "Finally, numpy does a lot of smart broadcasting of arrays. Broadcasting is what means that summing two arrays gives you a third array that has elements that are each the sum of the relevant elements in the two original arrays. Put another way, it's what causes `x + y = z` (for arrays x and y with the same shape) to result in an array z for which $z_{ij} = x_{ij} + y_{ij}$.\n",
    "\n",
    "Summing two arrays of the same shape is a pretty obvious example, but it also applies to cases that are *not* completely matched. For example, multiplication by a scalar is broadcast across all elements of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7d525c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.ones(shape=(3, 3))\n",
    "x * 3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e4646540",
   "metadata": {},
   "source": [
    "Similarly, numpy functions are broadcast across elements of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b9f50d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "828eab3b",
   "metadata": {},
   "source": [
    "This has been the slightest introduction to arrays. For a more comprehensive overview of arrays in Python, check out the excellent [\"NumPy: the absolute basics for beginners\"](https://numpy.org/doc/stable/user/absolute_beginners.html) on the **numpy** website. If you still want more delicious arrays, check out the [scipy lecture notes](https://scipy-lectures.org/intro/numpy/index.html)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d7b9f044",
   "metadata": {},
   "source": [
    "## Numeric linear algebra\n",
    "\n",
    "The transpose of an array `x` is given by `x.T`.\n",
    "\n",
    "Matrix multiplication is performed using the `@` operator. Here we perform $ M_{il} = \\sum_{k} x_{ik} * (x^T)_{kl}$, where $x^T$ is the transpose of $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0c6178",
   "metadata": {},
   "outputs": [],
   "source": [
    "x @ x.T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7b33fcc2",
   "metadata": {},
   "source": [
    "To multiply two arrays element wise, ie to do $ M_{ij} = x_{ij} y_{ij}$, it's the usual multiplication operator `*`.\n",
    "\n",
    "Inverting matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0e2478",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(9, size=(3, 3), dtype=\"int\")\n",
    "b = a @ np.linalg.inv(a)\n",
    "b"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ba258519",
   "metadata": {},
   "source": [
    "Computing the trace:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3907537a",
   "metadata": {},
   "outputs": [],
   "source": [
    "b.trace()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f5246ba0",
   "metadata": {},
   "source": [
    "Determinant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d3a6e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.det(a)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0c80fee0",
   "metadata": {},
   "source": [
    "Computing a Cholesky decomposition, i.e. finding lower triangular matrix $C$ such that $C C' = \\Sigma$ for $\\Sigma$ a 2-dimensional positive definite matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370af25f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Σ = np.array([[4, 1], [1, 3]])\n",
    "c = np.linalg.cholesky(Σ)\n",
    "\n",
    "c @ c.T - Σ"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9aeb0d0a",
   "metadata": {},
   "source": [
    "### Solving systems of linear equations\n",
    "\n",
    "Say we have a system of equations, $4x + 3y + 2z = 25$,\n",
    "$-2x + 2y + 3z = -10$, and $3x -5y + 2z = -4$. We can solve these three equations for the three unknowns, x, y, and z, using the `solve` method. First, remember that this equation can be written in matrix form as\n",
    "\n",
    "$$\n",
    "M\\cdot \\vec{x} = \\vec{c}\n",
    "$$\n",
    "\n",
    "We can solve this by multiplying by the matrix inverse of $M$:\n",
    "\n",
    "$$\n",
    "M^{-1} M \\cdot \\vec{x} = I \\cdot \\vec{x} = M^{-1} \\cdot \\vec{c}\n",
    "$$\n",
    "\n",
    "which could be called by running `x = la.inv(M).dot(c)`. There's a convenience function in **numpy** called solve that does the same thing: here it finds the real values of the vector $\\vec{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fcce7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[4, 3, 2], [-2, 2, 3], [3, -5, 2]])\n",
    "c = np.array([25, -10, -4])\n",
    "np.linalg.solve(M, c)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "dcdedda6",
   "metadata": {},
   "source": [
    "Finally, eigenvalues and eigenvectors can be found from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b82e345f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as la\n",
    "\n",
    "eigvals, eigvecs = la.eig(M)\n",
    "eigvals"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ad28b24c",
   "metadata": {},
   "source": [
    "## Interpolation\n",
    "\n",
    "This section draws on the **scipy** documentation. There are built-in **pandas** methods for interpolation in dataframes, but **scipy** also has a range of functions for this including for univariate data `interp1d`, multidimensional interpolation on a grid `interpn`, `griddata` for unstructured data. Let's see a simple example with interpolation between a regular grid of integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107bbe55",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "\n",
    "x = np.arange(0, 10)\n",
    "y = np.exp(-x / 3.0)\n",
    "f = interpolate.interp1d(x, y, kind=\"cubic\")\n",
    "# Create a finer grid to interpolation function f\n",
    "xnew = np.arange(0, 9, 0.1)\n",
    "ynew = f(xnew)\n",
    "plt.plot(x, y, \"o\", xnew, ynew, \"-\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "10c3cef0",
   "metadata": {},
   "source": [
    "What about unstructured data? Let's create a Cobb-Douglas function on a detailed grid but then only retain a random set of the established points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3f41ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import griddata\n",
    "\n",
    "\n",
    "def cobb_doug(x, y):\n",
    "    alpha = 0.8\n",
    "    return x ** (alpha) * y ** (alpha - 1)\n",
    "\n",
    "\n",
    "# Take some random points of the Cobb-Douglas function\n",
    "points = np.random.rand(1000, 2)\n",
    "values = cobb_doug(points[:, 0], points[:, 1])\n",
    "\n",
    "# Create a grid\n",
    "grid_x, grid_y = np.mgrid[0.01:1:200j, 0.01:1:200j]\n",
    "\n",
    "# Interpolate the points we have onto the grid\n",
    "interp_data = griddata(points, values, (grid_x, grid_y), method=\"cubic\")\n",
    "\n",
    "# Plot results\n",
    "fig, axes = plt.subplots(1, 2)\n",
    "# Plot function & scatter of random points\n",
    "axes[0].imshow(\n",
    "    cobb_doug(grid_x, grid_y).T,\n",
    "    extent=(0, 1, 0, 1),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"plasma_r\",\n",
    "    vmin=0,\n",
    "    vmax=1,\n",
    ")\n",
    "axes[0].plot(points[:, 0], points[:, 1], \"r.\", ms=1.2)\n",
    "axes[0].set_title(\"Original + points\")\n",
    "# Interpolation of random points\n",
    "axes[1].imshow(\n",
    "    interp_data.T, extent=(0, 1, 0, 1), origin=\"lower\", cmap=\"plasma_r\", vmin=0, vmax=1\n",
    ")\n",
    "axes[1].set_title(\"Cubic interpolation\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "09c7e3f7",
   "metadata": {},
   "source": [
    "## Optimisation\n",
    "\n",
    "**scipy** has functions for minimising scalar functions, minimising multivariate functions with complex surfaces, and root-finding. Let's see an example of finding the minimum of a scalar function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "309ca60c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return x**2 + 10 * np.sin(x) - 1.2\n",
    "\n",
    "\n",
    "result = optimize.minimize(f, x0=0)\n",
    "result"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "22c47af2",
   "metadata": {},
   "source": [
    "The result of the optimisation is in the 'x' attribute of `result`. Let's see this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05db5051",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-10, 10, 0.1)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, f(x))\n",
    "ax.scatter(result.x, f(result.x), s=150, color=\"k\")\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"f(x)\", rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b4ae23e1",
   "metadata": {},
   "source": [
    "In higher dimensions, the minimisation works in much the same way, with the same function `optimize.minimize`. There are a LOT of minimisation options that you can pass to the `method=` keyword; the default is intelligently chosen from BFGS, L-BFGS-B, or SLSQP, depending upon whether you supply constraints or bounds.\n",
    "\n",
    "Root finding, aka solving equations of the form $f(x)=0$, is also catered for by **scipy**, through `optimize.root`. It works in much the same way as `optimizer.minimize`.\n",
    "\n",
    "In both of these cases, be warned that multiple roots and multiple minima can be hard to detect, and you may need to carefully specify the bounds or the starting positions in order to find the root you're looking for. Also, both of these methods can accept the Jacobian of the function you're working with as an argument, which is likely to improve performance with some solvers.\n",
    "\n",
    "## Numerical Integration\n",
    "\n",
    "**scipy** provides routines to numerically evaluate integrals in `scipy.integrate`, which you can find the documentation for [here](https://docs.scipy.org/doc/scipy/reference/integrate.html). Let's see an example using the 'vanilla' integration method, `quad()`, to solve a known function between given (numerical) limits:\n",
    "\n",
    "$$\n",
    "\\displaystyle\\int_0^{\\pi} \\sin(x) d x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c075633f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "res, err = quad(np.sin, 0, np.pi)\n",
    "res"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f0c74920",
   "metadata": {},
   "source": [
    "What if we just have data samples? In that case, there are several routines that perform purely numerical integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c416809",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import simpson\n",
    "\n",
    "x = np.arange(0, 10)\n",
    "f_of_x = np.arange(0, 10)\n",
    "\n",
    "simpson(f_of_x, x) - 9**2 / 2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e021e70a",
   "metadata": {},
   "source": [
    "Even with just 10 evenly spaced points, the composite Simpson’s rule integration given by `simps` is able to accurately find the answer as $\\left( x^2/2\\right) |_{0}^{9}$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3ec3b8a2",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "In this chapter, you should have:\n",
    "\n",
    "- ✅ seen how to code numerical mathematics, including linear algebra and optimisation."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "formats": "md:myst"
  },
  "kernelspec": {
   "display_name": "codeforecon",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.16"
  },
  "vscode": {
   "interpreter": {
    "hash": "c4570b151692b3082981c89d172815ada9960dee4eb0bedb37dc10c95601d3bd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}