# Chapter 5
# Numerical Integration and Differentiation

In many computational economic applications, one must compute the definite integral
of a real-valued function f with respect to a "weighting" function w over an interval
$I$ of $R^n$:

$$\int_I f(x)w(x) dx$$


The weighting function may be the identity, $w = 1$, in which case the integral represents
the area under the function f. In other applications, w may be the probability
density of a random variable $\tilde X$
, in which case the integral represents the expectation
of $f( \tilde X)$ when $I$ repesents the whole support of $\tilde X$.






```python

```

In this chapter, we discuss three classes of numerical integration or numerical
quadrature methods<sup>1</sup>. All methods approximate the integral with a weighted sum of
function values:

$$\int_I f(x) w(x)dx \approx \sum_{i=0}^{n} w_i f(x_i)\thinspace .$$

<sup>1</sup>Quadrature is a historical mathematical term that means calculating area. 

The methods differ only in how the *quadrature weights* $wi$ and the *quadrature nodes*
$xi$ are chosen.

**Newton-Cotes** methods approximate the integrand f between nodes
using low order polynomials, and sum the integrals of the polynomials to estimate
the integral of f. Newton-Cotes methods are easy to implement, but are not particularly
eÆcient for computing the integral of a smooth function.

**Gaussian quadrature**
methods choose the nodes and weights to satisfy moment matching conditions, and
are more powerful than Newton-Cotes methods if the integrand is smooth.

**Monte Carlo and quasi-Monte Carlo integration** methods use "random" or "equidistributed"
nodes, and are simple to implement and are especially useful if the integration domain
is of high dimension or irregularly shaped.


In this chapter, we also present an overview of how to compute *finite difference*
approximations for the derivatives of a real-valued function. 

As we have seen in
previous chapters, it is often desirable to compute derivatives numerically because
analytic derivative expressions are difficult or impossible to derive, or expensive to
evaluate. 

Finite difference methods can also be used to solve differential equations,
which arise frequently in dynamic economic models, especially models formulated in
continuous time. In this chapter, we introduce numerical methods for differential
equations and illustrate their application to *initial value problems*.


```python
# https://github.com/QuantEcon/QuantEcon.py/blob/488b7b3b9117cfd9bfc71c187efc87c39fc5b459/quantecon/quad.py
"""
Filename: quad.py
Authors: Chase Coleman, Spencer Lyon
Date: 2014-07-01
Defining various quadrature routines.
Based on the quadrature routines found in the CompEcon toolbox by
Miranda and Fackler.
References
----------
Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
and Finance, MIT Press, 2002.
"""

from __future__ import division

import math
import numpy as np
import scipy.linalg as la
from scipy.special import gammaln
import sympy as sym
#from .ce_util import ckron, gridmake


from functools import reduce
```


```python
def ckron(*arrays):
    """
    Repeatedly applies the np.kron function to an arbitrary number of
    input arrays
    Parameters
    ----------
    *arrays : tuple/list of np.ndarray
    Returns
    -------
    out : np.ndarray
        The result of repeated kronecker products
    Notes
    -----
    Based of original function `ckron` in CompEcon toolbox by Miranda
    and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    return reduce(np.kron, arrays)

def gridmake(*arrays):
    """
    TODO: finish this docstring
    Notes
    -----
    Based of original function ``gridmake`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
    and Finance, MIT Press, 2002.
    """
    if all([i.ndim == 1 for i in arrays]):
        d = len(arrays)
        if d == 2:
            out = _gridmake2(*arrays)
        else:
            out = _gridmake2(arrays[0], arrays[1])
            for arr in arrays[2:]:
                out = _gridmake2(out, arr)

        return out
    else:
        raise NotImplementedError("Come back here")

        
def _gridmake2(x1, x2):
    """
    TODO: finish this docstring
    Notes
    -----
    Based of original function ``gridmake2`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
    and Finance, MIT Press, 2002.
    """
    if x1.ndim == 1 and x2.ndim == 1:
        return np.column_stack([np.tile(x1, x2.shape[0]),
                               np.repeat(x2, x1.shape[0])])
    elif x1.ndim > 1 and x2.ndim == 1:
        first = np.tile(x1, (x2.shape[0], 1))
        second = np.repeat(x2, x1.shape[0])
        return np.column_stack([first, second])
    else:
        raise NotImplementedError("Come back here")        

def _qnwtrap1(n, a, b):
    """
    Compute univariate trapezoid rule quadrature nodes and weights
    Parameters
    ----------
    n : int
        The number of nodes
    a : int
        The lower endpoint
    b : int
        The upper endpoint
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        An n element array of nodes
    nodes : np.ndarray(dtype=float)
        An n element array of weights
    Notes
    -----
    Based of original function ``qnwtrap1`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    if n < 1:
        raise ValueError("n must be at least one")

    nodes = np.linspace(a, b, n)
    dx = nodes[1] - nodes[0]

    weights = dx * np.ones(n)
    weights[0] *= 0.5
    weights[-1] *= 0.5

    return nodes, weights
        

def _qnwsimp1(n, a, b):
    """
    Compute univariate Simpson quadrature nodes and weights
    Parameters
    ----------
    n : int
        The number of nodes
    a : int
        The lower endpoint
    b : int
        The upper endpoint
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        An n element array of nodes
    nodes : np.ndarray(dtype=float)
        An n element array of weights
    Notes
    -----
    Based of original function ``qnwsimp1`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    if n % 2 == 0:
        print("WARNING qnwsimp: n must be an odd integer. Increasing by 1")
        n += 1

    nodes = np.linspace(a, b, n)
    dx = nodes[1] - nodes[0]
    weights = np.tile([2.0, 4.0], (n + 1) // 2)
    weights = weights[:n]
    weights[0] = weights[-1] = 1
    weights = (dx / 3.0) * weights

    return nodes, weights    
    
        
def _qnwlege1(n, a, b):
    """
    Compute univariate Guass-Legendre quadrature nodes and weights
    Parameters
    ----------
    n : int
        The number of nodes
    a : int
        The lower endpoint
    b : int
        The upper endpoint
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        An n element array of nodes
    nodes : np.ndarray(dtype=float)
        An n element array of weights
    Notes
    -----
    Based of original function ``qnwlege1`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    # import ipdb; ipdb.set_trace()
    maxit = 100
    m = np.fix((n + 1) / 2.0).astype(int)
    xm = 0.5 * (b + a)
    xl = 0.5 * (b - a)
    nodes = np.zeros(n)

    weights = nodes.copy()
    i = np.arange(m, dtype='int')

    z = np.cos(np.pi * ((i + 1.0) - 0.25) / (n + 0.5))

    for its in range(maxit):
        p1 = 1.0
        p2 = 0.0
        for j in range(1, n+1):
            p3 = p2
            p2 = p1
            p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j

        pp = n * (z * p1 - p2)/(z * z - 1.0)
        z1 = z.copy()
        z = z1 - p1/pp
        if all(np.abs(z - z1) < 1e-14):
            break

    if its == maxit - 1:
        raise ValueError("Maximum iterations in _qnwlege1")

    nodes[i] = xm - xl * z
    nodes[- i - 1] = xm + xl * z

    weights[i] = 2 * xl / ((1 - z * z) * pp * pp)
    weights[- i - 1] = weights[i]

    return nodes, weights





        
def _make_multidim_func(one_d_func, n, *args):
    """
    A helper function to cut down on code repetition. Almost all of the
    code in qnwcheb, qnwlege, qnwsimp, qnwtrap is just dealing
    various forms of input arguments and then shelling out to the
    corresponding 1d version of the function.
    This routine does all the argument checking and passes things
    through the appropriate 1d function before using a tensor product
    to combine weights and nodes.
    Parameters
    ----------
    one_d_func : function
        The 1d function to be called along each dimension
    n : int or array_like(float)
        A length-d iterable of the number of nodes in each dimension
    args :
        These are the arguments to various qnw____ functions.  For the
        majority of the functions this is just a and b, but some differ.
    Returns
    -------
    func : function
        The multi-dimensional version of the parameter ``one_d_func``
    """
    args = list(args)
    n = np.asarray(n)
    args = list(map(np.asarray, args))

    if all([x.size == 1 for x in [n] + args]):
        return one_d_func(n, *args)

    d = n.size

    for i in range(len(args)):
        if args[i].size == 1:
            args[i] = np.repeat(args[i], d)

    nodes = []
    weights = []

    for i in range(d):
        ai = [x[i] for x in args]
        _1d = one_d_func(n[i], *ai)
        nodes.append(_1d[0])
        weights.append(_1d[1])

    weights = ckron(*weights[::-1])  # reverse ordered tensor product

    nodes = gridmake(*nodes)
    return nodes, weights    
```

## 5.1 Newton-Cotes Methods



Newton-Cotes quadrature methods are designed to approximate the integral of a realvalued
function $f$ defined on a bounded interval $[a; b]$ of the real line. Newton-Cotes
methods approximate the integrand $f$ between nodes using *low order polynomials*,
and sum the integrals of the polynomials to form an estimate the integral of f. 

Two
Newton-Cotes rules are widely used in practice: the **trapezoid rule and Simpson's
rule**. Both rules are very easy to implement and are typically adequate for computing
the area under a continuous function.


```python

```

The trapezoid rule partitions the interval [a; b] into subintervals of equal length, approximates
f over each subinterval using linear interpolants, and then sums the areas under the
linear segments. The trapezoid rule draws its name from the fact that the area under f is
approximated by a series of trapezoids.


```python

```

where $x_i = a + (i-1)h$, with $h$ (called the step size) equal to   $ h=(b − a) / (n-1)$. The $w_i$ are called weights.

$$\int _{a}^{b}f(x)\,dx\approx \sum _{{i=1}}^{{n-1}}w_{i}\,f(x_{i}).$$


where $w_1 = w_n = h/2$ and $w_i = h$, otherwise.



```latex
%%latex
\begin{align}
\int_a^b f(x)\,dx &= \int_{x_0}^{x_1} f(x) dx + \int_{x_1}^{x_2} f(x) dx + \ldots + \int_{x_{n-1}}^{x_n} f(x) dx,     \nonumber \\ 
                  &\approx h \frac{f(x_0) + f(x_1)}{2} +
		  h \frac{f(x_1) + f(x_2)}{2} + \ldots + \nonumber \\ 
		  &\quad h \frac{f(x_{n-1}) + f(x_n)}{2} 
\end{align}
```


\begin{align}
\int_a^b f(x)\,dx &= \int_{x_0}^{x_1} f(x) dx + \int_{x_1}^{x_2} f(x) dx + \ldots + \int_{x_{n-1}}^{x_n} f(x) dx,     \nonumber \\ 
                  &\approx h \frac{f(x_0) + f(x_1)}{2} +
		  h \frac{f(x_1) + f(x_2)}{2} + \ldots + \nonumber \\ 
		  &\quad h \frac{f(x_{n-1}) + f(x_n)}{2} 
\end{align}




$$\int_a^b f(x)\,dx \approx  
\frac{h}{2}\left[f(x_0) + 2 f(x_1) + 2 f(x_2) + \ldots + 2 f(x_{n-1}) + f(x_n)\right]              
$$




$$
\int_a^b f(x)\,dx \approx h \left[\frac{1}{2}f(x_0) + \sum_{i=1}^{n-1}f(x_i) + \frac{1}{2}f(x_n) \right] \thinspace .
$$

For example, when $n = 2$


$${\frac  {b-a}{2}}(f_{0}+f_{1})$$


```python
def trapezoidal(f, a, b, n):
    h = float(b-a)/n
    result = 0.5*f(a) + 0.5*f(b)
    for i in range(1, n):
        result += f(a + i*h)
    result *= h
    return result
```


```python

```

The trapezoid rule is simple and robust. It is said to be first order exact because, if
not for rounding error, it will exactly compute the integral of any first order polynomial,
that is, a line. In general, if the integrand f is smooth, the trapezoid rule will yield an
approximation error that is $O(h^2)$, that is, the error shrinks quadratically with the width of
the subintervals.


```python
def qnwtrap(n, a, b):
    """
    Computes multivariate trapezoid rule quadrature nodes and weights.
    Parameters
    ----------
    n : int or array_like(float)
        A length-d iterable of the number of nodes in each dimension
    a : scalar or array_like(float)
        A length-d iterable of lower endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    b : scalar or array_like(float)
        A length-d iterable of upper endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        Quadrature nodes
    weights : np.ndarray(dtype=float)
        Weights for quadrature nodes
    Notes
    -----
    Based of original function ``qnwtrap`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    return _make_multidim_func(_qnwtrap1, n, a, b)
```


```python

```


```python

```

Simpson's rule is based on piece-wise quadratic, rather than piece-wise linear, approximations
to the integrand $f$.

$$\int_a^b f(x)dx \approx \sum_{i=0}^{n-1} w_if(x_i)\thinspace .$$

More formally, let $x_i = a + (i - 1)h$ for $i = 1; 2; ... ; n$, where
$ h=(b − a) / (n-1)$ and $n$ is odd. The nodes $x_i$ divide the interval $[a; b]$ into an even number
$n - 1$ of subintervals of equal length $h$.


```python

```


```python
def qnwsimp(n, a, b):
    """
    Computes multivariate Simpson quadrature nodes and weights.
    Parameters
    ----------
    n : int or array_like(float)
        A length-d iterable of the number of nodes in each dimension
    a : scalar or array_like(float)
        A length-d iterable of lower endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    b : scalar or array_like(float)
        A length-d iterable of upper endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        Quadrature nodes
    weights : np.ndarray(dtype=float)
        Weights for quadrature nodes
    Notes
    -----
    Based of original function ``qnwsimp`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    return _make_multidim_func(_qnwsimp1, n, a, b)
```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python
#https://github.com/birocoles/Disciplina-metodos-computacionais/blob/master/Content/newton-cotes.ipynb




# http://nbviewer.jupyter.org/github/sbustamante/ComputationalMethods/blob/master/material/numerical-calculus.ipynb#Numerical-Integration
```


```python

```


```python

```


```python

```

## 5.2 Gaussian Quadrature
Gaussian quadrature rules are constructed with respect to specific weighting functions.
Specifically, for a weighting function $w$ defined on an interval $I \in R$ of the real
line, and for a given order of approximation n, the quadrature nodes $x_1; x_2; ... ; x_n$
and quadrature weights $w_1; w_2; ...; w_n$ are chosen so as to satisfy the $2n$ "momentmatching"
conditions:


```python

```


```python

```


```python

```

## 5.3 Monte Carlo Integration
Monte Carlo integration methods are motivated by the Strong Law of Large Numbers.
One version of the Law states that if $x_1; x_2; ... ; x_n$ are independent realizations of a
random variable $\tilde X$
and $f$ is a continuous function, then




```python

```


```python

```

## 5.4 Quasi-Monte Carlo Integration
Although Monte-Carlo integration methods originated using insights from probability
theory, recent extensions have severed that connection and, in the process, demonstrated
ways in which the methods can be improved. Monte-Carlo methods rely on
sequences $x_i$ with the property that


```python

```


```python

```

## 5.5 An Integration Toolbox
The Matlab toolbox accompanying the textbook includes four functions for computing
numerical integrals for general functions. Each takes three inputs, n, a, and b and
generates appropriate nodes and weights. The functions `qnwtrap` and `qnwsimp` implement the Newton-Cotes trapezoid and Simpson's rule methods, `qnwlege` implements
Gauss-Legendre quadrature and `qnwequi` generates nodes and weights associated with
either equidistributed or pseudo-random sequences. The calling syntax is the same
for each and is illustrated with below with `qnwtrap`.





```python




    
        







def qnwlege(n, a, b):
    """
    Computes multivariate Guass-Legendre  quadrature nodes and weights.
    Parameters
    ----------
    n : int or array_like(float)
        A length-d iterable of the number of nodes in each dimension
    a : scalar or array_like(float)
        A length-d iterable of lower endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    b : scalar or array_like(float)
        A length-d iterable of upper endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        Quadrature nodes
    weights : np.ndarray(dtype=float)
        Weights for quadrature nodes
    Notes
    -----
    Based of original function ``qnwlege`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    return _make_multidim_func(_qnwlege1, n, a, b)



def qnwequi(n, a, b, kind="N", equidist_pp=None):
    """
    Generates equidistributed sequences with property that averages
    value of integrable function evaluated over the sequence converges
    to the integral as n goes to infinity.
    Parameters
    ----------
    n : int
        Number of sequence points
    a : scalar or array_like(float)
        A length-d iterable of lower endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    b : scalar or array_like(float)
        A length-d iterable of upper endpoints. If a scalar is given,
        that constant is repeated d times, where d is the number of
        dimensions
    kind : string, optional(default="N")
        One of the following:
        - N - Neiderreiter (default)
        - W - Weyl
        - H - Haber
        - R - pseudo Random
    equidist_pp : array_like, optional(default=None)
        TODO: I don't know what this does
    Returns
    -------
    nodes : np.ndarray(dtype=float)
        Quadrature nodes
    weights : np.ndarray(dtype=float)
        Weights for quadrature nodes
    Notes
    -----
    Based of original function ``qnwequi`` in CompEcon toolbox by
    Miranda and Fackler
    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.
    """
    if equidist_pp is None:
        equidist_pp = np.sqrt(np.array(list(sym.primerange(0, 7920))))

    n, a, b = list(map(np.atleast_1d, list(map(np.asarray, [n, a, b]))))

    d = max(list(map(len, [n, a, b])))
    n = np.prod(n)

    if a.size == 1:
        a = np.repeat(a, d)

    if b.size == 1:
        b = np.repeat(b, d)

    i = np.arange(1, n + 1)

    if kind.upper() == "N":  # Neiderreiter
        j = 2.0 ** (np.arange(1, d+1) / (d+1))
        nodes = np.outer(i, j)
        nodes = (nodes - np.fix(nodes)).squeeze()
    elif kind.upper() == "W":  # Weyl
        j = equidist_pp[:d]
        nodes = np.outer(i, j)
        nodes = (nodes - np.fix(nodes)).squeeze()
    elif kind.upper() == "H":  # Haber
        j = equidist_pp[:d]
        nodes = np.outer(i * (i+1) / 2, j)
        nodes = (nodes - np.fix(nodes)).squeeze()
    elif kind.upper() == "R":  # pseudo-random
        nodes = np.random.rand(n, d).squeeze()
    else:
        raise ValueError("Unknown sequence requested")

    # compute nodes and weights
    r = b - a
    nodes = a + nodes * r
    weights = (np.prod(r) / n) * np.ones(n)

    return nodes, weights



```


```python

```


```python

```


```python

```


```python

```

All of the quadrature functions will use tensor products to generate nodes and
weights for integration over an arbitrary bounded interval [a; b] in higher dimensional
spaces.


```python

```


```python

```

## 5.6 Numerical Differentiation
The most natural way to approximate a derivative is to replace it with a finite difference.
The definition of a derivative,





```python

```


```python

```


```python

```

## 5.7 Initial Value Problems
Differential equations pose the problem of inferring a function given information
about its derivatives and additional "boundary" conditions. Differential equations
may characterized as either ordinary differential equations (ODEs), whose solutions
are functions of a single argument, and partial differential equations (PDEs), whose
solutions are functions of multiple arguments. Both ODEs and PDEs may be solved
numerically using finite difference methods.

From a numerical point of view the distinction between ODEs and PDEs is less
important than the distinction between initial value problems (IVPs), which can be
solved in a recursive or evolutionary fashion, and boundary value problems (BVPs),
which require the entire solution to be computed simultaneously because the solution
at one point (in time and/or space) depends on the solution everywhere else. For
ODEs, the solution of an IVP is known at some point and the solution near this
point can then be (approximately) determined. This, in turn, allows the solution at
still other points to be approximated and so forth. BVPs, on the other hand, require
simultaneous solution of the differential equation and the boundary conditions. We
take up the solution of IVPs in this section, but defer discussion of BVPs until the
next chapter (page 164).


```python

```


```python

```

There are numerous other approaches and refinements to solving initial value problems.
Brie
y, these include so-called multi-step algorithms which utilize information
from previous steps to determine the current step direction (Runge-Kutta are singlestep
methods). Also, any method can adapt the step size to the current behavior
of the system by monitoring the truncation error, reducing (increasing) the step size
if this error is unacceptably large (small). Adaptive schemes are important if one
requires a given level of accuracy.


```python

```

### Example: Commercial Fishery


```python

```


```python

```


```python

```


```python

```

b3c74322e8d98bf2e89a951e6e2b6a87d1c06c2d
ipynb
Jupyter Notebook
Chapter05.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
2018-05-09T08:17:44.000Z
2021-12-26T07:02:17.000Z
Chapter05.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
Chapter05.ipynb
lnsongxf/Applied_Computational_Economics_and_Finance
f14661bfbfa711d49539bda290d4be5a25087185
MIT
2017-12-15T13:39:35.000Z
2021-05-15T15:06:02.000Z
Qwen/Qwen-72B
1. YES
2. YES