

# <span style="color:#2c061f"> Exercise 11 </span>  

<br>

## <span style="color:#374045"> Introduction to Programming and Numerical Analysis </span>


####  <span style="color:#d89216"> <br> Sebastian Honoré </span>

## Plan for today
<br>

1. Welcome
2. Linear algebra in python 
3. Solving equations analytically 
4. Problemset 6

**Remember peer feedback deadline is Sunday 24th by 23.59**

# 2.1 Linear algebra in python 

We can use `SciPy's` module `linalg` to perform a range of operations within the realm of linear algebra: 

- Basic stuff: determinant, invert, norm
- Matrix decompositions (LU, Cholesky etc.)
- Solve a system of equations
- Find eigenvalues

# 2.2 Solve a system of equations
Lets take the module for a spin. We want to solve:

$$Ax=b$$


```python
import numpy as np 
from scipy import linalg
np.random.seed(2021)
A = np.random.uniform(size=(5,5))
b = np.random.uniform(size=5)
print(f'Matrix A: \n{A} \n\n Matrix b:\n{b}')
```

    Matrix A: 
    [[0.60597828 0.73336936 0.13894716 0.31267308 0.99724328]
     [0.12816238 0.17899311 0.75292543 0.66216051 0.78431013]
     [0.0968944  0.05857129 0.96239599 0.61655744 0.08662996]
     [0.56127236 0.61652471 0.96384302 0.57430429 0.37116085]
     [0.45214524 0.20185025 0.56930512 0.19509597 0.58370402]] 
    
     Matrix b:
    [0.47631347 0.5178144  0.82309863 0.73222503 0.06905627]


# 2.3 Solve a system of equations

There are several ways we can solve this system of equations:


```python
# Using a regular solver ()
x = linalg.solve(A,b)
x
```




    array([ 3.17983125, -1.88260751, -2.05057824,  4.3774516 , -1.15692713])




```python
# Using LU decomposition 
#-> decompose A into upper and lower triangular matrix and solve through substitution
LU, piv = linalg.lu_factor(A)
x = linalg.lu_solve((LU,piv),b)
x
```




    array([ 3.17983125, -1.88260751, -2.05057824,  4.3774516 , -1.15692713])



# 2.4 When to use solve and lu_factorization?

Combined linalg.lu_factor and linalg.lu_solve is equivalent to linalg.solve. But the algorithmic complexity is different! 

- linalg.solve is $O(n^3)$, but linalg.lu_solve is $O(n^2)$
- So when only the b matrix changes and the A matrix is constant it is more efficient to use linalg.lu_solve. 
- If the matrix is high dimensional linalg.lu_solve is fastest. 

# 2.4 Why care? 

Being able to solve a system of linear equations probably won't earn much street credit from your friends. However, it may come in handy. Remember the good ol' OLS equation? In case you forgot: 

$$\hat{\beta} = (X'X)^{-1}X'y$$

In the 2020 exam you had to code the estimator from scratch and solve it using the functionalities of linalg. How?


```python
# Solving OLS 
def beta_hat(x1,x2,y,con=True):
    """ 
    Estimates the vector of coefficients for two independent variables
    
    Args:
    x1 (array): Array of independent variable values
    x2 (array): Array of independent variable values
    y (array): Array of dependent variable values
    con (bool): Include array of intercepts. Defaults to True. 

    Returns:
    beta_hat (float): vector of estimated coefficients
    """
    
    if con == True:
        con = np.ones((x1.shape)) #if con=True create a N x 1 vector, for the intercepts to be estimated
        X = np.stack((con,x1,x2),axis=1) #Stack vectors into the matrix X, as shown above
    else:
        X = np.stack((x1,x2),axis=1) #Stack vectors into the matrix X, as shown above
    
    beta_hat = linalg.inv(X.T@X)@X.T@y #The matrix of coeffients is computed given the definition
    return beta_hat
```

# 3.1 Solving equations analytically 
Formulas and models are part of the life as an economist. `Sympy` enables us to translate these formulas into python code that we can work with - super cool! Lets see it in action:

Consider a utility function from a standard OLG model. Economic agents lives two periods (young/old) and obtains utility from consumption in both periods:

$$U_{i,t} = u_i(c_{i,1}) + \frac{1}{1+\rho}u_i(C_{i,2})$$

# 3.2 Implementation in Sympy
For the sake of simplicity let's assume log-utility:


```python
import sympy as sm
c1,c2 = sm.symbols("C_i1"), sm.symbols("C_i2")
rho = sm.symbols("rho")

#utility
uc1 = sm.ln(c1)
uc2 = sm.ln(c2)
U = uc1+1/(1+rho)*uc2
U
```




$\displaystyle \log{\left(C_{i1} \right)} + \frac{\log{\left(C_{i2} \right)}}{\rho + 1}$



# 3.3 Sympy functionalities
We can do a lot of stuff from here. Say we need the derivative of $U$ with respect to $C_{i,2}$


```python
sm.diff(U,c2)
```




$\displaystyle \frac{1}{C_{i2} \left(\rho + 1\right)}$



**Differentation has never been easier**

# 3.4 Turn your formulas into functions
A cool feature of Sympy is that you can turn your formulas into python functions. This will likely be useful in your model project:

- Solve model analytically with Sympy
- Convert solution into python function e.g. the law-of-motion in the OLG model
- Find steady state using an optimizer

How is it done?

# 3.5 Functions in Sympy
 
Use the `Lambdify` method which takes a function and an iterable as argument. In our case the function is utility and the iterable is consumption: 


```python
util = sm.lambdify((c1,c2,rho),U)



print(f'Utility for C1=10, C2=5 and rho=0.1: {util(10,5,0.1):.2f}')
```

    Utility for C1=10, C2=5 and rho=0.1: 3.77


# Problemset 6

Your turn to work with `linalg` and `SymPy`.

d8edd2c37532495da4fa34c3555c3bfe73d32641
ipynb
Jupyter Notebook
class 11/exercise 11.ipynb
s-honore/NumEcon_slides22
3ad4027c721dd3c3d4fe5c62cf3bb22a1dffeb08
MIT
2022-02-08T10:44:33.000Z
2022-02-08T10:44:33.000Z
class 11/exercise 11.ipynb
s-honore/NumEcon_slides22
3ad4027c721dd3c3d4fe5c62cf3bb22a1dffeb08
MIT
class 11/exercise 11.ipynb
s-honore/NumEcon_slides22
3ad4027c721dd3c3d4fe5c62cf3bb22a1dffeb08
MIT
2022-02-09T10:52:28.000Z
2022-02-09T10:52:28.000Z
Qwen/Qwen-72B
1. YES
2. YES