# Matrix and compartive statistics review

The following notebook is a review of matrices and comparative statistics with examples in python.

The equations and examples are from the following book I highly recommend using to brush up on mathamtics commonly used in economics coursework:
- Dowling, E. T. (2012). Introduction to mathematical economics. McGraw-Hill.
    - [Amazon link](https://www.amazon.com/Schaums-Introduction-Mathematical-Economics-Outlines/dp/0071762515/ref=sr_1_7?dchild=1&keywords=mathematics+economics&qid=1593200726&sr=8-7)

# Table of contents
- [1. Matrix basics](#1.-Matrix-basics)
    
- [2. Special determinants](#2.-Special-determinants)

- [3. Comparative statistics](#3.-Comparative-statistics)



# 1. Matrix basics


```python
import numpy as np
np.random.seed(1)
```

## 1.1 Scalar multiplication


```python
A = np.random.randint(20, size=(2,2))
A
```




    array([[ 5, 11],
           [12,  8]])




```python
A*3
```




    array([[15, 33],
           [36, 24]])



## 1.2 Matrix addition


```python
A = np.random.randint(20, size=(2,2))
B = np.random.randint(20, size=(2,2))
A+B
```




    array([[ 9, 27],
           [ 6, 27]])



## 1.3 Matrix mulitiplication


```python
A = np.random.randint(20, size=(2,2))
B = np.random.randint(20, size=(2,2))
```


```python
A@B
```




    array([[178, 256],
           [228, 288]])



### 1.4 Identity matrix and empty


```python
np.eye(3)
```




    array([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.]])




```python
np.zeros(shape=(2,2))
```




    array([[0., 0.],
           [0., 0.]])



## 1.5 Matrix inversion


```python
A = np.random.randint(20, size=(2,2))
A
```




    array([[14, 18],
           [ 4,  9]])




```python
#Determinate
round(np.linalg.det(A))
```




    54.0




```python
#Invert matrix
np.linalg.inv(A).round(3)
```




    array([[ 0.167, -0.333],
           [-0.074,  0.259]])



# 2. Special determinants

## 2.1 Jacobian


```python
import sympy as sy
x1, x2 = sy.symbols('x1 x2', integer=True)
y1 = 5*x1+3*x2
y2 = 25*x1**2+30*x1*x2+9*x2**2
```


```python
y1
```




$\displaystyle 5 x_{1} + 3 x_{2}$




```python
y2
```




$\displaystyle 25 x_{1}^{2} + 30 x_{1} x_{2} + 9 x_{2}^{2}$




```python
independent_variables = [x1, x2]
functions = [y1, y2]
Jacobian = sy.Matrix(np.zeros(shape=(len(functions), len(independent_variables))))
count = 0
for funcs in functions:
    for iv in independent_variables:
        Jacobian[count] = sy.diff(funcs, iv)
        count+=1
Jacobian
```




$\displaystyle \left[\begin{matrix}5 & 3\\50 x_{1} + 30 x_{2} & 30 x_{1} + 18 x_{2}\end{matrix}\right]$




```python
if Jacobian.det()==0:
    print("Functional dependence")
```

    Functional dependence


## 2.2 Hessian

### 2.2.1 Sympy example: Definition


```python
from sympy import Function, hessian
from sympy.abc import x, y
f = Function('f')(x, y)
hessian(f, (x, y))
```




$\displaystyle \left[\begin{matrix}\frac{\partial^{2}}{\partial x^{2}} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)}\\\frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y^{2}} f{\left(x,y \right)}\end{matrix}\right]$



### 2.2.1 Sympy example: From 2.1


```python
count = 0
Hessian = Jacobian.copy()
for _ in range(0,2):
    for iv in independent_variables: #Reverses list order for hessian
        Hessian[count] = sy.diff(Hessian[count],iv)
        count+=1
Hessian
```




$\displaystyle \left[\begin{matrix}0 & 0\\50 & 18\end{matrix}\right]$




```python
H1 = Hessian[0]
H2 = Hessian.det()
if H1>0 and H2>0:
    print('Positive definite')
    print('Minimum point')
if H1<0 and H2>0:
    print('Negative definite')
    print('Max point')
```

## 2.3 Discriminant
- Tests for positive or negative definiteness of quadratic equations


```python
x, y = sy.symbols('x y', integer=True)
z = 2*x**2 + 5*x*y+8*y**2
z
```




$\displaystyle 2 x^{2} + 5 x y + 8 y^{2}$




```python
Discrim = sy.Matrix([[2,(5/2)], [(5/2),8]])
Discrim 
```




$\displaystyle \left[\begin{matrix}2 & 2.5\\2.5 & 8\end{matrix}\right]$




```python
d1 = Discrim[0]
d2 = Discrim.det() #deter
if d1 and d2>0:
    print("Positive definite")
elif d1<0 and d2>0:
    print('Negative definite')
```

    Positive definite


## 2.4 Higher order hessian


```python
x1, x2, x3 = sy.symbols('x1 x2 x3', integer=True)
z = -5*x1**2+10*x1+x1*x3-2*x2**2+4*x2+2*x2*x3-4*x3**2
z
```




$\displaystyle - 5 x_{1}^{2} + x_{1} x_{3} + 10 x_{1} - 2 x_{2}^{2} + 2 x_{2} x_{3} + 4 x_{2} - 4 x_{3}^{2}$




```python
focs = []
for idx,iv in enumerate([x1,x2,x3]):
    foc = sy.diff(z,iv)
    focs.append(foc)
for idx, foc in enumerate(focs):
    print("FOC %s:" %(idx+1), foc)
```

    FOC 1: -10*x1 + x3 + 10
    FOC 2: -4*x2 + 2*x3 + 4
    FOC 3: x1 + 2*x2 - 8*x3



```python
A,b = sy.linear_eq_to_matrix(focs, [x1, x2, x3])
```


```python
Hessian = A
Hessian
```




$\displaystyle \left[\begin{matrix}-10 & 0 & 1\\0 & -4 & 2\\1 & 2 & -8\end{matrix}\right]$




```python
H1 = Hessian[0]
H2 = Hessian[0:2,0:2].det()
H3 = Hessian.det()
if H1>0 and H2>0 and H3>0:
    print('Positive definite')
    print('Minimum point')
elif H1<0 and H2>0 and H3<0:
    print('Negative definite')
    print('Maximum point')
```

    Negative definite
    Maximum point


## 2.5 Bordered Hessian


```python
from sympy import Function, hessian, pprint
from sympy.abc import x, y
f = Function('f')(x, y)
constraint = Function('g')(x, y)
hessian(f, (x, y), [constraint])
```




$\displaystyle \left[\begin{matrix}0 & \frac{\partial}{\partial x} g{\left(x,y \right)} & \frac{\partial}{\partial y} g{\left(x,y \right)}\\\frac{\partial}{\partial x} g{\left(x,y \right)} & \frac{\partial^{2}}{\partial x^{2}} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)}\\\frac{\partial}{\partial y} g{\left(x,y \right)} & \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)} & \frac{\partial^{2}}{\partial y^{2}} f{\left(x,y \right)}\end{matrix}\right]$



## 2.6 Eigenvalues & Eigenvectors


```python
def eigen(matrix):
    trace = np.trace(A)
    det = round(np.linalg.det(A),0)
    eig_values = np.round((np.sort((trace+np.array([+1,-1])*np.sqrt(trace**2-(4*det)))/2)),1)
    solu1, solu2 = [eig_values[:][0], eig_values[:][1]]
    print("Original matrix: \n",matrix)
    print("Eigen-values:\n {}".format(eig_values))
    #Classification of matrix
    if solu1>0 and solu2>0:
        print('Pos definite')
    if solu1<0 and solu2<0:
        print('Neg definite')
    if (solu1==0 or solu2==0) and (solu1>0 and solu2>0):
        print('Pos semi-def')
    if (solu1==0 or solu2==0) and (solu1<0 and solu2<0):
        print('Pos semi-def')
    if (solu1<0 and solu2>0) or (solu1>0 and solu2<0):
        print('Indefinite')
```


```python
A = np.random.randint(20, size=(2,2))
eigen(A)
```

    Original matrix: 
     [[17  0]
     [13  9]]
    Eigen-values:
     [ 9. 17.]
    Pos definite


# 3. Comparative statistics

## 3.1 One endogenous variable

$$
Q_d = m-nP+kY\\
Q_s = a+bP
$$

### 3.1.1 Explicit function
$$P* =\frac{m-a+kY}{b+n}$$

### 3.1.2 Implicit function

$$\frac{dP^*}{dY}= - \frac{F_Y}{F_P}$$


```python
from sympy.abc import x,n,p,k,y,a,b,m
f = m-n*p+k*y-a-b*p
f
```




$\displaystyle - a - b p + k y + m - n p$




```python
-sy.diff(f,y)/sy.diff(f,p)
```




$\displaystyle - \frac{k}{- b - n}$



## 3.2 N-endogenous variables
- `Comparative statistics:` requires a unique equilibrium condition for each endogenous variable
- Measuring the effect of an exogenous variable on the endgenous variables involves taking the total derivative of each equilibrium conditions
    - w.r.t to the particular exogenous variable and solving for each of the partial derivatives

$$
F^1(y_1, y_2; x_1, x_2) = 0 \\
F^2(y_1, y_2; x_1, x_2) = 0
$$

#### Note: 
- #### Exogenous variables: $x_1$ and $x_2$ 
- #### Endogenous variables: $y_1$ and $y_2$ 

## 3.3 Comparative statistics for optimization problems
- Apply comparative statistics to the first order conditions to determine initial optimal values


```python
from sympy.abc import r, K, w, L, P, Q
Q = Function('Q')(K, L)
π = p*Q-r*K-w*L
π
```




$\displaystyle - K r - L w + p Q{\left(K,L \right)}$



### 3.3.1 F.O.C


```python
focs = []
for idx,iv in enumerate([K,L]):
    foc = sy.diff(π,iv)
    focs.append(foc)
focs[0]
```




$\displaystyle p \frac{\partial}{\partial K} Q{\left(K,L \right)} - r$



### 3.3.2 Jacobian

- For optimization of a system the Deteriminant of the Jacobian>0


```python
Jacobian = sy.Matrix([[π.diff(K,K), π.diff(K,L)],[π.diff(L,K), π.diff(L,L)]])
Jacobian
```




$\displaystyle \left[\begin{matrix}p \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} & p \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}\\p \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & p \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)}\end{matrix}\right]$




```python
B = []
for foc in focs:
    B.append(foc.diff(r))
B
```




    [-1, 0]



### 3.3.3 Find derivatives


```python
J = Jacobian.det()
```


```python
J1 = Jacobian.copy()
J1[0] =1
J1[2] =0
J1
```




$\displaystyle \left[\begin{matrix}1 & p \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}\\0 & p \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)}\end{matrix}\right]$




```python
J2 = Jacobian.copy()
J2[1] =1
J2[3] =0
J2
```




$\displaystyle \left[\begin{matrix}p \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} & 1\\p \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & 0\end{matrix}\right]$



#### 3.3.3.1 Find $\frac{\partial \bar{K}}{\partial r}$


```python
J1.det()/J
```




$\displaystyle \frac{p \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)}}{p^{2} \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} - p^{2} \left(\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}\right)^{2}}$



#### 3.3.3.2 Find $\frac{\partial \bar{L}}{\partial r}$


```python
J2.det()/J
```




$\displaystyle - \frac{p \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}}{p^{2} \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} - p^{2} \left(\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}\right)^{2}}$



## 3.4 Comparative statistics in constrained optimization
- Optimize comparative statistics with constraints


```python
from sympy.abc import r, K, w, L, P, Q, B
lamda = sy.symbols('lamda')
Q = Function('Q')(K, L)
π = Q+lamda*(B-r*K-w*L)
π 
```




$\displaystyle \lambda \left(B - K r - L w\right) + Q{\left(K,L \right)}$



### 3.4.1 F.O.C


```python
focs = []
for idx,iv in enumerate([K,L,lamda]):
    foc = sy.diff(π,iv)
    focs.append(foc)
```


```python
focs[0]
```




$\displaystyle - \lambda r + \frac{\partial}{\partial K} Q{\left(K,L \right)}$




```python
focs[1]
```




$\displaystyle - \lambda w + \frac{\partial}{\partial L} Q{\left(K,L \right)}$




```python
focs[2]
```




$\displaystyle B - K r - L w$



### 3.4.2 Jacobian


```python
independent_variables = [K,L,lamda]
functions = focs
Jacobian = sy.Matrix(np.zeros(shape=(len(functions), len(independent_variables))))
count = 0
for funcs in functions:
    for iv in independent_variables:
        Jacobian[count] = sy.diff(funcs, iv)
        count+=1
Jacobian
```




$\displaystyle \left[\begin{matrix}\frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} & \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & - r\\\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} & - w\\- r & - w & 0\end{matrix}\right]$




```python
J_deter = Jacobian.det()
```

### 3.4.3 Find derivatives



```python
def deriv_convert(matrix, col=0):
    deriv_col = iter([0, 0, -1])
    derivs = Jacobian.copy()
    for i in range(1,len(Jacobian)+1):
        if i%3==col:
            derivs[i-1] = next(deriv_col)
    return derivs
```

#### 3.3.4.1 Find $\frac{\partial \bar{K}}{\partial B}$


```python
k_b = deriv_convert(derivs, col=1)
k_b
```




$\displaystyle \left[\begin{matrix}0 & \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & - r\\0 & \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} & - w\\-1 & - w & 0\end{matrix}\right]$




```python
k_b.det()/J_deter
```




$\displaystyle \frac{- r \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} + w \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}}{- r^{2} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} + 2 r w \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} - w^{2} \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)}}$



#### 3.3.4.2 Find $\frac{\partial \bar{L}}{\partial B}$


```python
L_b = deriv_convert(derivs, col=2)
L_b
```




$\displaystyle \left[\begin{matrix}\frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} & 0 & - r\\\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & 0 & - w\\- r & -1 & 0\end{matrix}\right]$




```python
L_b.det()/J_deter
```




$\displaystyle \frac{r \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} - w \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)}}{- r^{2} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} + 2 r w \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} - w^{2} \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)}}$



#### 3.3.4.3 Find $\frac{\partial \bar{\lambda}}{\partial B}$


```python
lamb_b = deriv_convert(derivs, col=0)
lamb_b 
```




$\displaystyle \left[\begin{matrix}\frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} & \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & 0\\\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} & \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} & 0\\- r & - w & -1\end{matrix}\right]$




```python
lamb_b.det()/J_deter
```




$\displaystyle \frac{- \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} + \left(\frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)}\right)^{2}}{- r^{2} \frac{\partial^{2}}{\partial L^{2}} Q{\left(K,L \right)} + 2 r w \frac{\partial^{2}}{\partial L\partial K} Q{\left(K,L \right)} - w^{2} \frac{\partial^{2}}{\partial K^{2}} Q{\left(K,L \right)}}$



## 3.5 Envelope theorem

- `Envelope theorem:` Measures the effect of a change in exogenous variables on the optimal value of the objective function
    - This can be achieved by simply taking the derivative of the Lagrangian function w.r.t the desired exogenous variable and evaluating the derivative at the values of the optimal solution
    


```python
B, x, px, y, py, lamb= sy.symbols('B x px y py lamda', integer=True)

u = Function('u')(x, y)
constraint = lamb*(B-px*x-py*y)
U = u+constraint
print('Budget constraint:')
U
```

    Budget constraint:





$\displaystyle \lambda \left(B - px x - py y\right) + u{\left(x,y \right)}$




```python
focs = []
for idx,iv in enumerate([px, py, B]):
    foc = sy.diff(U,iv)
    focs.append(foc)
```

- $\lambda$: Marginal uility of money
    - The extra utility derived from a change in income
    
- The first and second order conditions are negative:
    - $\uparrow \lambda \rightarrow$  negative impact on the quantity of good consumed


```python
focs[0]
```




$\displaystyle - \lambda x$




```python
focs[1]
```




$\displaystyle - \lambda y$




```python
focs[2]
```




$\displaystyle \lambda$



## 3.6 Concave programming
- Optimize comparative statistics with inequality constraints
- Assume that functions are concave


```python
B, x, px, y, py, lamb= sy.symbols('B x px y py lamda', integer=True)

u = Function('u')(x, y)
constraint = lamb*(B-px*x-py*y)
U = u+constraint
print('Budget constraint:')
U
```

    Budget constraint:





$\displaystyle \lambda \left(B - px x - py y\right) + u{\left(x,y \right)}$




```python
#1.A
sy.diff(U, x)
```




$\displaystyle - \lambda px + \frac{d}{d x} u{\left(x,y \right)}$




```python
# 1.B
sy.diff(U, y)
```




$\displaystyle - \lambda py + \frac{d}{d y} u{\left(x,y \right)}$




```python
#2.A
sy.diff(U,lamb)
```




$\displaystyle B - px x - py y$



dc4b1b4917cf1a7defba70b6b2b51fa45fa5cc31
ipynb
Jupyter Notebook
Math_reviews/Matrix_ComparativeStats_Review.ipynb
wenddymacro/DSGE_models-1
0e3bba5ab9c80ac954a78fd9d8d8a0930e47b2f7
MIT
2020-07-14T05:56:35.000Z
2022-01-04T10:14:41.000Z
Math_reviews/Matrix_ComparativeStats_Review.ipynb
corybaird/DSGE_models
ccb4ecdda8944a532423687b28191c8bd3c61612
MIT
Math_reviews/Matrix_ComparativeStats_Review.ipynb
corybaird/DSGE_models
ccb4ecdda8944a532423687b28191c8bd3c61612
MIT
2021-01-17T16:46:35.000Z
2022-03-26T12:28:29.000Z
Qwen/Qwen-72B
1. YES
2. YES