```python
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 3: NumPy

NumPy is a powerful Python module for scientific computing. Among other things, NumPy defines an N-dimensional array object that is especially convenient to use for plotting functions and for simulating and storing time series data. NumPy also defines many useful mathematical functions like, for example, the sine, cosine, and exponential functions and has excellent functions for probability and statistics including random number generators, and many cumulative density functions and probability density functions.

## Importing NumPy

The standard way to import NumPy so that the namespace is `np`. This is for the sake of brevity.


```python
import numpy as np
```

## NumPy arrays

A NumPy `ndarray` is a homogeneous multidimensional array. Here, *homogeneous* means that all of the elements of the array have the same type. An `nadrray` is a table of numbers (like a matrix but with possibly more dimensions) indexed by a tuple of positive integers. The dimensions of NumPy arrays are called axes and the number of axes is called the rank. For this course, we will work almost exclusively with 1-dimensional arrays that are effectively vectors. Occasionally, we might run into a 2-dimensional array.


### Basics

The most straightforward way to create a NumPy array is to call the `array()` function which takes as an argument a `list`. For example:


```python
# Create a variable called a1 equal to a numpy array containing the numbers 1 through 5
a1 = np.array([1,2,3,4,5])
print(a1)

# Find the type of a1
print(type(a1))

# Find the shape of a1
print(np.shape(a1))

# Use ndim to find the rank or number of dimensions of a1
print(np.ndim(a1))
```

    [1 2 3 4 5]
    <class 'numpy.ndarray'>
    (5,)
    1



```python
# Create a variable called a2 equal to a 2-dimensionl numpy array containing the numbers 1 through 4
a2 = np.array([[1,2],[3,4]])
print(a2)

# Use ndim to find the rank or number of dimensions of a2
print(np.ndim(a2))

# Find the shape of a2
print(np.shape(a2))
```

    [[1 2]
     [3 4]]
    2
    (2, 2)



```python
# Create a variable called c an empty numpy array
a3 = np.array([])
print(a3)

# Use ndim to find the rank or number of dimensions of a3
print(np.ndim(a3))

# Find the shape of a3
print(np.shape(a3))
```

    []
    1
    (0,)


### Special functions for creating arrays

Numpy has several built-in functions that can assist you in creating certain types of arrays: `arange()`, `zeros()`, and `ones()`. Of these, `arrange()` is probably the most useful because it allows you a create an array of numbers by specifying the initial value in the array, the maximum value in the array, and a step size between elements. `arrange()` has three arguments: `start`, `stop`, and `step`:

    arange([start,] stop[, step,])
  
The `stop` argument is required. The default for `start` is 0 and the default for `step` is 1. Note that the values in the created array will stop one increment *below* stop. That is, if `arrange()` is called with `stop` equal to 9 and `step` equal to 0.5, then the last value in the returned array will be 8.5.


```python
# Create a variable called b that is equal to a numpy array containing the numbers 1 through 5
b = np.arange(1,6,1)
print(b)
```

    [1 2 3 4 5]



```python
# Create a variable called c that is equal to a numpy array containing the numbers 0 through 10
c = np.arange(11)
print(c)
```

    [ 0  1  2  3  4  5  6  7  8  9 10]


The `zeros()` and `ones()` take as arguments the desired shape of the array to be returned and fill that array with either zeros or ones.


```python
# Construct a 1x5 array of zeros
print(np.zeros(5))
```

    [0. 0. 0. 0. 0.]



```python
# Construct a 2x2 array of ones
print(np.zeros([2,2]))
```

    [[0. 0.]
     [0. 0.]]


### Math with NumPy arrays

A nice aspect of NumPy arrays is that they are optimized for mathematical operations. The following standard Python arithemtic operators `+`, `-`, `*`, `/`, and `**` operate element-wise on NumPy arrays as the following examples indicate.


```python
# Define three 1-dimensional arrays. PROVIDED
A = np.array([2,4,6])
B = np.array([3,2,1])
C = np.array([-1,3,2,-4])
```


```python
# Multiply A by a constant
print(3*A)
```

    [ 6 12 18]



```python
# Exponentiate A
print(A**2)
```

    [ 4 16 36]



```python
# Add  A and B together
print(A+B)
```

    [5 6 7]



```python
# Exponentiate A with B
print(A**B)
```

    [ 8 16  6]



```python
# Add A and C together
print(A+C)
```

The error in the preceding example arises because addition is element-wise and `A` and `C` don't have the same shape.


```python
# Compute the sine of the values in A
print(np.sin(A))
```

    [ 0.90929743 -0.7568025  -0.2794155 ]


### Iterating through Numpy arrays

NumPy arrays are iterable objects just like lists, strings, tuples, and dictionaries which means that you can use `for` loops to iterate through the elements of them.


```python
# Use a for loop with a NumPy array to print the numbers 0 through 4
for x in np.arange(5):
    print(x)
```

    0
    1
    2
    3
    4


### Example: Basel problem

One of my favorite math equations is:

\begin{align}
\sum_{n=1}^{\infty} \frac{1}{n^2} & = \frac{\pi^2}{6}
\end{align}

We can use an iteration through a NumPy array to approximate the lefthand-side and verify the validity of the expression.


```python
# Set N equal to the number of terms to sum
N = 1000

# Initialize a variable called summation equal to 0
summation = 0

# loop over the numbers 1 through N
for n in np.arange(1,N+1):
    summation = summation + 1/n**2

# Print the approximation and the exact solution
print('approx:',summation)
print('exact: ',np.pi**2/6)
```

    approx: 1.6439345666815615
    exact:  1.6449340668482264


## Random numbers and statistics

NumPy has many useful routines for generating draws from probability distributions. Random number generation is useful for computing dynamic simulations of stochastic economies. NumPy also has some routines for computing basic summary statistics (min, max, mean, median, variance) for data in an array.

For more advanced applications involving probability distributions, use SciPy's `scipy.stats` module (https://docs.scipy.org/doc/scipy-0.16.1/reference/stats.html). And for estimating statisticsl models (e.g., OLS models) use StatsModels (https://www.statsmodels.org/stable/index.html).


```python
# Use the help function to view documentation on the np.random.uniform function
help(np.random.uniform)
```

    Help on built-in function uniform:
    
    uniform(...) method of mtrand.RandomState instance
        uniform(low=0.0, high=1.0, size=None)
        
        Draw samples from a uniform distribution.
        
        Samples are uniformly distributed over the half-open interval
        ``[low, high)`` (includes low, but excludes high).  In other words,
        any value within the given interval is equally likely to be drawn
        by `uniform`.
        
        Parameters
        ----------
        low : float or array_like of floats, optional
            Lower boundary of the output interval.  All values generated will be
            greater than or equal to low.  The default value is 0.
        high : float or array_like of floats
            Upper boundary of the output interval.  All values generated will be
            less than high.  The default value is 1.0.
        size : int or tuple of ints, optional
            Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
            ``m * n * k`` samples are drawn.  If size is ``None`` (default),
            a single value is returned if ``low`` and ``high`` are both scalars.
            Otherwise, ``np.broadcast(low, high).size`` samples are drawn.
        
        Returns
        -------
        out : ndarray or scalar
            Drawn samples from the parameterized uniform distribution.
        
        See Also
        --------
        randint : Discrete uniform distribution, yielding integers.
        random_integers : Discrete uniform distribution over the closed
                          interval ``[low, high]``.
        random_sample : Floats uniformly distributed over ``[0, 1)``.
        random : Alias for `random_sample`.
        rand : Convenience function that accepts dimensions as input, e.g.,
               ``rand(2,2)`` would generate a 2-by-2 array of floats,
               uniformly distributed over ``[0, 1)``.
        
        Notes
        -----
        The probability density function of the uniform distribution is
        
        .. math:: p(x) = \frac{1}{b - a}
        
        anywhere within the interval ``[a, b)``, and zero elsewhere.
        
        When ``high`` == ``low``, values of ``low`` will be returned.
        If ``high`` < ``low``, the results are officially undefined
        and may eventually raise an error, i.e. do not rely on this
        function to behave when passed arguments satisfying that
        inequality condition.
        
        Examples
        --------
        Draw samples from the distribution:
        
        >>> s = np.random.uniform(-1,0,1000)
        
        All values are within the given interval:
        
        >>> np.all(s >= -1)
        True
        >>> np.all(s < 0)
        True
        
        Display the histogram of the samples, along with the
        probability density function:
        
        >>> import matplotlib.pyplot as plt
        >>> count, bins, ignored = plt.hist(s, 15, density=True)
        >>> plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
        >>> plt.show()
    



```python
# Print a random draw from the uniform(0,1) probability distribution.
print(np.random.uniform())
```

    0.39395193886380386


In the previous example, everyone will obtain different output for `np.random.uniform()`. However, we can set the *seed* for NumPy's PRNG (pseudorandom number generator) with the function `np.random.seed()`. Setting the seed before generating a random sequence can ensure replicability of results.


```python
# Set the seed for the NumPy random number generator to 271828
np.random.seed(271828)

# Print a random draw from the uniform(0,1) probability distribution.
print(np.random.uniform())
```

    0.9662022955536614


### Example: Random Draws from $\text{uniform}(-1,1)$ Distribution

Create a sample of 200 draws from the $\text{uniform}(-1,1)$ distribution. Compute some summary statistics. Plot the sample.


```python
# Set the seed for the NumPy random number generator to 314159
np.random.seed(314159)

# Create a variable called 'uniform_data' containing 200 random draw from the uniform(-1,1) probability distribution.
uniform_data = np.random.uniform(size=200,low=-1,high=1) 

# Print the variable data
print(uniform_data)
```

    [ 0.63584662  0.10209259 -0.16044929 -0.80261629  0.6220415   0.93471281
     -0.80358661  0.60372074  0.20980423  0.16953762 -0.01382494 -0.04650081
     -0.47206371 -0.38079105  0.81590802 -0.63003033 -0.62352117 -0.82321921
      0.74842975  0.28481629  0.70565033  0.0661206   0.75880411 -0.41319904
     -0.14432861  0.38336681 -0.94102104 -0.73408651 -0.29515564  0.3641906
      0.9278488  -0.09180626  0.02434442 -0.06863382 -0.73856441 -0.58805703
     -0.53966153 -0.33441377  0.3820414  -0.43905549  0.65183529 -0.66164868
     -0.72066187  0.55183355 -0.52219358  0.26076713  0.89209106 -0.47310432
     -0.50685854  0.66302636  0.44291651 -0.26846445  0.81962684  0.9394919
      0.68700008  0.09848171  0.25145527  0.00982847  0.74560833  0.13846091
     -0.6741501  -0.78951496  0.67829627  0.77302013  0.92888613  0.01447009
     -0.33037939  0.45926088  0.93300648 -0.23495154  0.19028247 -0.71479702
      0.71114841  0.07243895 -0.632449   -0.69894595 -0.65026494  0.9580402
     -0.26946754  0.85194828  0.38758249  0.64178299 -0.6738458   0.83502118
      0.63629364  0.02859558  0.68737679 -0.21023369 -0.85377648 -0.84879093
      0.81242546 -0.94086188 -0.10959221 -0.01291769 -0.95679783  0.85599001
     -0.17968129  0.99989443  0.23891269 -0.66466892  0.9739976   0.95522212
      0.69345454  0.42483686  0.74119329 -0.25768943 -0.24289615  0.3512096
     -0.01420467 -0.81386935 -0.75655242  0.67474155  0.81256821  0.09794052
     -0.35859784  0.11619177  0.52145372 -0.26214142  0.35253046  0.55032724
      0.18885249  0.78057581  0.83120461 -0.45686194  0.8972115  -0.96391123
      0.89429982 -0.81937678 -0.52383     0.91979452 -0.62225414  0.8141226
      0.50637822 -0.6921192   0.77288675  0.27131884 -0.90608287  0.76175962
      0.90954684 -0.13271693 -0.81301378  0.60804403  0.92988145  0.39675914
      0.37276877 -0.57952684 -0.2294779  -0.86317973 -0.16198036  0.28673408
     -0.57565177 -0.94692889 -0.63431133 -0.49885546 -0.98798956  0.64392303
     -0.22363103  0.35045596  0.77442238  0.9363158  -0.90516324 -0.55774438
     -0.57081079 -0.13020764 -0.53775636 -0.22456364  0.24785222  0.90810479
      0.29336168 -0.38997457 -0.10495809 -0.13725783  0.44436708  0.2566703
      0.12444098 -0.7196481   0.65897494  0.92976743  0.31022742  0.17614419
     -0.96249917  0.88456632  0.39453716 -0.75154362 -0.02893635 -0.73586317
     -0.80650354  0.02100266  0.46247197 -0.88487504 -0.18950471 -0.70548356
     -0.63094273  0.45872944  0.70059834  0.67150071  0.51254239  0.48583715
     -0.52763227 -0.68046815]



```python
# Print the mean of the values in variable uniform_data
print('mean:  ',np.mean(uniform_data))

# Print the median of the values in variable uniform_data
print('median:',np.median(uniform_data))

# Print the standard deviation of the values in variable uniform_data
print('min:   ',np.std(uniform_data))

# Print the maximum of the values in variable uniform_data
print('max:   ',np.max(uniform_data))

# Print the index value of the arg max of the values in variable uniform_data
print('argmax:',np.argmax(uniform_data))
```

    mean:   0.047187009866624184
    median: 0.06927977272929209
    min:    0.6090136251058506
    max:    0.9998944319955057
    argmax: 97



```python
# Create a variable called index_max equal to the argmax of data. Print the corresponding element of uniform_data
index_max = np.argmax(uniform_data)
print('data at argmax:',uniform_data[index_max])
```

    data at argmax: 0.9998944319955057



```python
# Plot random sample
plt.plot(uniform_data)
plt.title('200 uniform(-1,1) draws')
plt.grid()
```

### Example: Random Draws from $\text{normal}\left(0,\frac{1}{3}\right)$ Distribution

Create a sample of 200 draws from a normal distribution with mean 0 and variance $\frac{1}{3}$. Note that the variance of the $\text{uniform}(-1,1)$ distribution from the previous example has the same variance See: https://en.wikipedia.org/wiki/Uniform_distribution_(continuous).

Compute some summary statistics. Plot the sample.


```python
# Use the help function to view documentation on the np.random.normal function
help(np.random.normal)
```

    Help on built-in function normal:
    
    normal(...) method of mtrand.RandomState instance
        normal(loc=0.0, scale=1.0, size=None)
        
        Draw random samples from a normal (Gaussian) distribution.
        
        The probability density function of the normal distribution, first
        derived by De Moivre and 200 years later by both Gauss and Laplace
        independently [2]_, is often called the bell curve because of
        its characteristic shape (see the example below).
        
        The normal distributions occurs often in nature.  For example, it
        describes the commonly occurring distribution of samples influenced
        by a large number of tiny, random disturbances, each with its own
        unique distribution [2]_.
        
        Parameters
        ----------
        loc : float or array_like of floats
            Mean ("centre") of the distribution.
        scale : float or array_like of floats
            Standard deviation (spread or "width") of the distribution.
        size : int or tuple of ints, optional
            Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
            ``m * n * k`` samples are drawn.  If size is ``None`` (default),
            a single value is returned if ``loc`` and ``scale`` are both scalars.
            Otherwise, ``np.broadcast(loc, scale).size`` samples are drawn.
        
        Returns
        -------
        out : ndarray or scalar
            Drawn samples from the parameterized normal distribution.
        
        See Also
        --------
        scipy.stats.norm : probability density function, distribution or
            cumulative density function, etc.
        
        Notes
        -----
        The probability density for the Gaussian distribution is
        
        .. math:: p(x) = \frac{1}{\sqrt{ 2 \pi \sigma^2 }}
                         e^{ - \frac{ (x - \mu)^2 } {2 \sigma^2} },
        
        where :math:`\mu` is the mean and :math:`\sigma` the standard
        deviation. The square of the standard deviation, :math:`\sigma^2`,
        is called the variance.
        
        The function has its peak at the mean, and its "spread" increases with
        the standard deviation (the function reaches 0.607 times its maximum at
        :math:`x + \sigma` and :math:`x - \sigma` [2]_).  This implies that
        `numpy.random.normal` is more likely to return samples lying close to
        the mean, rather than those far away.
        
        References
        ----------
        .. [1] Wikipedia, "Normal distribution",
               http://en.wikipedia.org/wiki/Normal_distribution
        .. [2] P. R. Peebles Jr., "Central Limit Theorem" in "Probability,
               Random Variables and Random Signal Principles", 4th ed., 2001,
               pp. 51, 51, 125.
        
        Examples
        --------
        Draw samples from the distribution:
        
        >>> mu, sigma = 0, 0.1 # mean and standard deviation
        >>> s = np.random.normal(mu, sigma, 1000)
        
        Verify the mean and the variance:
        
        >>> abs(mu - np.mean(s)) < 0.01
        True
        
        >>> abs(sigma - np.std(s, ddof=1)) < 0.01
        True
        
        Display the histogram of the samples, along with
        the probability density function:
        
        >>> import matplotlib.pyplot as plt
        >>> count, bins, ignored = plt.hist(s, 30, density=True)
        >>> plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
        ...                np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
        ...          linewidth=2, color='r')
        >>> plt.show()
    



```python
# Set the seed for the NumPy random number generator to 314159
np.random.seed(314159)

# Create a variable called 'normal_data' containing 200 random draw from the normal(0,1/3) probability distribution.
normal_data = np.random.normal(size=200,loc=0,scale=1/np.sqrt(3)) 

# Print the variable normal_data
print(normal_data)
```

    [ 0.12143549  0.75631678 -0.50674044 -0.10130139  0.83074966  1.0280597
     -1.92533743 -0.57241295 -0.51265814 -0.63553831  0.19357178  0.50866078
      0.06320691  0.67455494 -0.21112296  0.38770897  1.02091162 -0.38434928
      0.7808367  -0.63282346 -0.03008466  0.3040535  -1.76142292  0.62477673
     -0.17253536 -0.21669408 -0.40996677 -0.66158548 -0.6408805   0.55765817
     -0.22355639  0.22024066  0.21849396 -0.28534015  0.37852699 -0.75801105
      0.39010377 -0.29821956 -0.48548378  0.80095813  0.09902537  0.69079261
      0.0529736   1.35529688  0.11087969  0.5970842   0.00488096  0.3133261
      0.70744534 -0.50891633 -0.05539998  0.21999659 -0.61272051  0.1631092
      0.06779908  0.66559788 -0.20814988 -0.18834673  0.36935728 -0.1168261
      0.53048069  0.32036534  0.03481789  0.77474907 -0.19405353  0.63447439
     -0.20067492 -1.70250295 -0.64075609  0.23031732  0.27426255  0.44767446
     -0.18670597  0.53702323  0.87605589 -0.60587922 -0.52375627 -0.00914125
      0.0618483   0.51312738  0.35156255 -1.08501292 -0.38057956  0.75705177
      0.63411274  0.40620206  0.52555331  0.12715235 -0.12791087  0.23271823
     -0.36527786  0.26724985  0.17083956  0.48665855 -0.04839994  0.33169856
     -0.59265491  0.38121314 -0.37510574 -0.09972254  1.05956336 -0.59856316
     -0.33056808 -0.42032833  0.91195135 -0.58192938 -0.18789336 -0.82369635
     -0.32696855 -0.78298257  0.27393143  0.07476506 -0.78159695  0.58796295
     -1.21538751 -0.92938051  0.47173415  0.81670191 -0.63784809  0.11029619
      0.57882996  1.01944278 -0.41393313  0.2173021  -0.6381914  -0.0250956
      0.01392882 -0.53486761 -0.04042405  0.02112727 -0.62492897 -0.1678664
      0.33839992 -0.46543986  0.13840744  0.14440495  0.46850309  0.49425552
     -0.35292786 -0.27365884 -0.60058536 -0.24254545  0.68385251  0.29115578
      0.50685277 -0.4930715  -0.62612325  0.09506943  1.5985902  -0.44056757
     -1.07103108 -0.15234541 -0.37974221  0.02521733 -0.15139174  0.18952504
      0.45473891  0.04437202  0.46953245  0.40186251  0.57882478  0.03737229
     -0.51040583  0.4320896   0.13216004 -0.13552132  0.96610447  1.10346492
     -0.36091273 -0.03749187  0.06744852  0.06116802  0.35073528 -0.19900517
     -0.63715201  0.52139514  0.20782661 -0.018552    0.43281211 -0.55847675
     -0.30212557 -0.32978655 -0.35717671  0.93606909  0.59562641 -0.4201014
      0.12613396  0.208413   -0.48309414  0.05591969 -0.66945551 -0.39336486
      0.0247874  -0.49181131  0.63436971  0.91645032 -0.50152753 -0.0021372
      0.05158724 -0.14778507]



```python
# Print the mean of the values in variable normal_data
print('mean:  ',np.mean(normal_data))

# Print the median of variable the values in normal_data
print('median:',np.median(normal_data))

# Print the standard deviation of the values in variable normal_data
print('min:   ',np.std(normal_data))

# Print the maximum of the values in variable normal_data
print('max:   ',np.max(normal_data))

# Print the index value of the arg max of the values in variable normal_data
print('argmax:',np.argmax(normal_data))
```

    mean:   0.020031561013250983
    median: 0.040872154927040286
    min:    0.5552816540529858
    max:    1.5985901984025754
    argmax: 148



```python
# Create a variable called 'index_max' equal to the argmax of data. Print corresponding element of normal_data
index_max = np.argmax(normal_data)
print('data at argmax:',normal_data[index_max])
```

    data at argmax: 1.5985901984025754



```python
# Plot the random sample
plt.plot(normal_data)
plt.title('200 normal(0,1/3) draws')
plt.grid()
```


```python
# Plot the uniform and normal samples together on the same axes
plt.plot(normal_data,label='normal(0,1/3)')
plt.plot(uniform_data,label='uniform(-1,1)')
plt.legend(loc='lower right')
plt.title('Two random samples')
plt.grid()
```

8bd7faf7b35c2086f71d63fceb59ee9ba7dc0a62
ipynb
Jupyter Notebook
Lecture Notebooks/Econ126_Class_03.ipynb
pmezap/computational-macroeconomics
b703f46176bb16e712badf752784f8a7b996cdb1
MIT
2020-02-29T06:09:03.000Z
2022-03-25T13:14:13.000Z
Lecture Notebooks/Econ126_Class_03.ipynb
letsgoexploring/computational-macroeconomics
b703f46176bb16e712badf752784f8a7b996cdb1
MIT
Lecture Notebooks/Econ126_Class_03.ipynb
letsgoexploring/computational-macroeconomics
b703f46176bb16e712badf752784f8a7b996cdb1
MIT
2019-09-24T07:48:49.000Z
2022-03-22T21:36:30.000Z
Qwen/Qwen-72B
1. YES
2. YES