6.4. Applied Linear Algebra#

# import numpy to prepare for code below
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
# activate plot theme

6.4.1. Vectors and Matrices#

6.4.1.1. Vectors#

A (N-element) vector is \( N \) numbers stored together.

We typically write a vector as \( x = \begin{bmatrix} x_1 \\ x_2 \\ \dots \\ x_N \end{bmatrix} \).

In numpy terms, a vector is a 1-dimensional array.

In a previous lecture, we saw some types of operations that can be done on vectors, such as

x = np.array([1, 2, 3])
y = np.array([4, 5, 6])

Element-wise operations: Let \( z = x ? y \) for some operation \( ? \), one of the standard binary operations (\( +, -, \times, \div \)). Then we can write \( z = \begin{bmatrix} x_1 ? y_1 & x_2 ? y_2 \end{bmatrix} \). Element-wise operations require that \( x \) and \( y \) have the same size.

print("Element-wise Addition", x + y)
print("Element-wise Subtraction", x - y)
print("Element-wise Multiplication", x * y)
print("Element-wise Division", x / y)
Element-wise Addition [5 7 9]
Element-wise Subtraction [-3 -3 -3]
Element-wise Multiplication [ 4 10 18]
Element-wise Division [0.25 0.4  0.5 ]

Scalar operations: Let \( w = a ? x \) for some operation \( ? \), one of the standard binary operations (\( +, -, \times, \div \)). Then we can write \( w = \begin{bmatrix} a ? x_1 & a ? x_2 \end{bmatrix} \).

print("Scalar Addition", 3 + x)
print("Scalar Subtraction", 3 - x)
print("Scalar Multiplication", 3 * x)
print("Scalar Division", 3 / x)
Scalar Addition [4 5 6]
Scalar Subtraction [2 1 0]
Scalar Multiplication [3 6 9]
Scalar Division [3.  1.5 1. ]

Another operation very frequently used in data science is the dot product.

The dot between \( x \) and \( y \) is written \( x \cdot y \) and is equal to \( \sum_{i=1}^N x_i y_i \).

We can use @ to denote dot products (and matrix multiplication which we’ll see soon!).

print("Dot product with @", x @ y)
Dot product with @ 32

6.4.1.2. Matrices#

An \( N \times M \) matrix can be thought of as a collection of M N-element vectors stacked side-by-side as columns.

We write a matrix as

\[\begin{split} \begin{bmatrix} x_{11} & x_{12} & \dots & x_{1M} \\ x_{21} & \dots & \dots & x_{2M} \\ \vdots & \vdots & \vdots & \vdots \\ x_{N1} & x_{N2} & \dots & x_{NM} \end{bmatrix} \end{split}\]

In numpy terms, a matrix is a 2-dimensional array.

We can create a matrix by passing a list of lists to the np.array function.

x = np.array([[1, 2, 3], [4, 5, 6]])
y = np.ones((2, 3))
z = np.array([[1, 2], [3, 4], [5, 6]])

We can perform element-wise and scalar operations as we did with vectors. In fact, we can do these two operations on arrays of any dimension.

print("Element-wise Addition\n", x + y)
print("Element-wise Subtraction\n", x - y)
print("Element-wise Multiplication\n", x * y)
print("Element-wise Division\n", x / y)

print("Scalar Addition\n", 3 + x)
print("Scalar Subtraction\n", 3 - x)
print("Scalar Multiplication\n", 3 * x)
print("Scalar Division\n", 3 / x)
Element-wise Addition
 [[2. 3. 4.]
 [5. 6. 7.]]
Element-wise Subtraction
 [[0. 1. 2.]
 [3. 4. 5.]]
Element-wise Multiplication
 [[1. 2. 3.]
 [4. 5. 6.]]
Element-wise Division
 [[1. 2. 3.]
 [4. 5. 6.]]
Scalar Addition
 [[4 5 6]
 [7 8 9]]
Scalar Subtraction
 [[ 2  1  0]
 [-1 -2 -3]]
Scalar Multiplication
 [[ 3  6  9]
 [12 15 18]]
Scalar Division
 [[3.   1.5  1.  ]
 [0.75 0.6  0.5 ]]

Similar to how we combine vectors with a dot product, matrices can do what we’ll call matrix multiplication.

Matrix multiplication is effectively a generalization of dot products.

Matrix multiplication: Let \( v = x \cdot y \) then we can write \( v_{ij} = \sum_{k=1}^N x_{ik} y_{kj} \) where \( x_{ij} \) is notation that denotes the element found in the ith row and jth column of the matrix \( x \).

The image below from Wikipedia, by Bilou, shows how matrix multiplication simplifies to a series of dot products:

image

After looking at the math and image above, you might have realized that matrix multiplication requires very specific matrix shapes!

For two matrices \( x, y \) to be multiplied, \( x \) must have the same number of columns as \( y \) has rows.

Formally, we require that for some integer numbers, \( M, N, \) and \( K \) that if \( x \) is \( N \times M \) then \( y \) must be \( M \times K \).

If we think of a vector as a \( 1 \times M \) or \( M \times 1 \) matrix, we can even do matrix multiplication between a matrix and a vector!

Let’s see some examples of this.

x1 = np.reshape(np.arange(6), (3, 2))
x2 = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])
x3 = np.array([[2, 5, 2], [1, 2, 1]])
x4 = np.ones((2, 3))

y1 = np.array([1, 2, 3])
y2 = np.array([0.5, 0.5])

Numpy allows us to do matrix multiplication in three ways.

print("Using @ for two matrices")
print(x1 @ x4)
Using @ for two matrices
[[1. 1. 1.]
 [5. 5. 5.]
 [9. 9. 9.]]
print("Using @ for vec and mat")
print(y1 @ x1)
Using @ for vec and mat
[16 22]

6.4.2. Transpose#

A matrix transpose is an operation that flips all elements of a matrix along the diagonal.

More formally, the \( (i, j) \) element of \( x \) becomes the \( (j, i) \) element of \( x^T \).

In particular, let \( x \) be given by

\[\begin{split} x = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end{bmatrix} \end{split}\]

then \( x \) transpose, written as \( x' \), is given by

\[\begin{split} x = \begin{bmatrix} 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 9 \\ \end{bmatrix} \end{split}\]

In Python, we do this by

x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

print("x transpose is")
print(x.transpose())
x transpose is
[[1 4 7]
 [2 5 8]
 [3 6 9]]

6.4.3. Identity Matrix#

In linear algebra, one particular matrix acts very similarly to how 1 behaves for scalar numbers.

This matrix is known as the identity matrix and is given by

\[\begin{split} I = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ 0 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 1 \end{bmatrix} \end{split}\]

As seen above, it has 1s on the diagonal and 0s everywhere else.

When we multiply any matrix or vector by the identity matrix, we get the original matrix or vector back!

Let’s see some examples.

I = np.eye(3)
x = np.reshape(np.arange(9), (3, 3))
y = np.array([1, 2, 3])

print("I @ x", "\n", I @ x)
print("x @ I", "\n", x @ I)
print("I @ y", "\n", I @ y)
print("y @ I", "\n", y @ I)
I @ x 
 [[0. 1. 2.]
 [3. 4. 5.]
 [6. 7. 8.]]
x @ I 
 [[0. 1. 2.]
 [3. 4. 5.]
 [6. 7. 8.]]
I @ y 
 [1. 2. 3.]
y @ I 
 [1. 2. 3.]

6.4.4. Inverse#

If you recall, you learned in your primary education about solving equations for certain variables.

For example, you might have been given the equation

\[ 3x + 7 = 16 \]

and then asked to solve for \( x \).

You probably did this by subtracting 7 and then dividing by 3.

Now let’s write an equation that contains matrices and vectors.

\[\begin{split} \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 3 \\ 4 \end{bmatrix} \end{split}\]

How would we solve for \( x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \)?

Unfortunately, there is no “matrix divide” operation that does the opposite of matrix multiplication.

Instead, we first have to do what’s known as finding the inverse. We must multiply both sides by this inverse to solve.

Consider some matrix \( A \).

The inverse of \( A \), given by \( A^{-1} \), is a matrix such that \( A A^{-1} = I \) where \( I \) is our identity matrix.

Notice in our equation above, if we can find the inverse of \( \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix} \) then we can multiply both sides by the inverse to get

\[\begin{split} \begin{align*} \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}^{-1}\begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} &= \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}^{-1}\begin{bmatrix} 3 \\ 4 \end{bmatrix} \\ I \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} &= \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}^{-1} \begin{bmatrix} 3 \\ 4 \end{bmatrix} \\ \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} &= \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}^{-1} \begin{bmatrix} 3 \\ 4 \end{bmatrix} \end{align*} \end{split}\]

Computing the inverse requires that a matrix be square and satisfy some other conditions (non-singularity) that are beyond the scope of this lecture.

We also skip the exact details of how this inverse is computed, but, if you are interested, you can visit the QuantEcon Linear Algebra lecture for more details.

We demonstrate how to compute the inverse with numpy below.

# This is a square (N x N) non-singular matrix
A = np.array([[1, 2, 0], [3, 1, 0], [0, 1, 2]])

print("This is A inverse")

print(np.linalg.inv(A))

print("Check that A @ A inverse is I")
print(np.linalg.inv(A) @ A)
This is A inverse
[[-0.2  0.4  0. ]
 [ 0.6 -0.2  0. ]
 [-0.3  0.1  0.5]]
Check that A @ A inverse is I
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 2.77555756e-17  1.00000000e+00  0.00000000e+00]
 [-1.38777878e-17  0.00000000e+00  1.00000000e+00]]

6.4.5. Portfolios#

Vectors and inner products give us a convenient way to organize and calculate the value of portfolios of assets

6.4.5.1. Static Payoffs#

As an example, consider a portfolio with 4 dollars of asset A, 2.5 dollars of asset B, and 8 dollars of asset C in date 0.

The assets have the following (net) returns 3%, 5%, and -1.10% in date 1.

First, calculate the value of the portfolio in the end of date 1

4.0 *(1+3/100) + 2.5 * (1+5/100) + 8 * (1-1.1/100)
14.832999999999998

Attention

Why are we dividing by 100 and adding 1?

Because returns are a number like 5%. So If you had 100 dollars you end up with 105 dollars.

Meaning you get your “pricipal” plus the return.

So the “1” here is makign the role of the pricipal, so if you return is 0% you get 100.

The division by 100 is just making sure that this is in percentage, so that your gross retruns are a number like 1.05 when your net returns are 5%.

We can make this more convenient and general by using arrays for accounting, and then sum then in a loop.

import numpy as np
x = np.array([4.0, 2.5, 8.0]) # portfolio positio
y = np.array([3.0, 5.0, -1.1]) # Asset net returns 

n = len(x)
p = 0.0
for i in range(n): # i.e. 0, 1, 2
    p = p + x[i] * (1+y[i]/100)

p
14.832999999999998

The above would have worked with x and y as list rather than np.array.

Note that the general pattern above is the sum.

\[ p = \sum_{i=0}^{n-1} x_i y_i = x \cdot y \]

This is an inner product as implemented by the @ function

x @ (1+y/100)
14.832999999999998

This approach allows us to simultaneously price different portfolios by stacking them in a matrix and using the dot product (the @ operator does the dot product for us!).

y = np.array([3.0, 5.0, -1.1]) # returns
x1 = np.array([4.0, 2.5, 8.0]) # portfolio 1
x2 = np.array([2.0, 1.5, 0.0]) # portfolio 2
X = np.array((x1, x2))

# calculate with inner products
p1 = X[0,:] @ (1+y/100)
p2 = X[1,:] @ (1+y/100)
print("Calculating separately")
print([p1, p2])

# or with a matrix multiplication
print("Calculating with matrices")
P = X @ (1+y/100)
print(P)
Calculating separately
[14.832999999999998, 3.6350000000000002]
Calculating with matrices
[14.833  3.635]

Insight: If you substitute dividends by returns and number of shares by portfolio weights, this same formulation allow us to easily compute the returns of a portfolio

\[ r = \sum_{i=0}^{n-1} w_i r_i = w \cdot r \]

How do we construct the weights?

The value of each position divided by the toal value of the portfolio

What is the value of each position? x

What is the value of the porfolio in date 0? sum(x)

port_value=np.sum(x)
w=x/port_value

w
array([0.27586207, 0.17241379, 0.55172414])

Thus the net returns of the portfolio are

w@y
2.296551724137931

So the final value of the portfolio can be written as

port_value*(1+w@y/100)
14.833000000000002