

# <span style="color:#2c061f"> Exercise 10 </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. Data project revisited
3. Linear algebra in python 
4. Solving equations analytically 
5. Problemset 6

# 2.1 Data project revisited

Overall some very nice projects! 

**Things to keep in mind**
- Restart the kernel and run all cells before hand-in
- Clear stories. Give a small introduction to the problem at hand and how you solve it. 
- [The DRY principle](https://thevaluable.dev/dry-principle-cost-benefit-example/)
- Documentation, documentation, documentation! 
- Thoughtful plotting

# 3.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

# 3.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]


# 3.3 Solve a system of equations

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


```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])




```python
# Regular solver
linalg.solve(A,b)
```




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



# 3.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. 

# 4.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})$$

# 4.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}$



# 4.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**

# 4.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?

# 4.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


# 5 Problemset 6

Your turn to work with `linalg` and `SymPy`. Please go to your channels. 

98d119f729fc6f48cab03a63e5c28ad5abb751c3
ipynb
Jupyter Notebook
Slides/exercise 10/exercise 10.ipynb
s-honore/NumEcon
e720a1eef6993150868d69a90791b8ff4054f6b3
MIT
2021-02-24T13:01:19.000Z
2021-02-24T13:01:19.000Z
Slides/exercise 10/exercise 10.ipynb
s-honore/NumEcon
e720a1eef6993150868d69a90791b8ff4054f6b3
MIT
Slides/exercise 10/exercise 10.ipynb
s-honore/NumEcon
e720a1eef6993150868d69a90791b8ff4054f6b3
MIT
2021-02-24T14:41:20.000Z
2021-02-24T14:41:20.000Z
Qwen/Qwen-72B
1. YES
2. YES