Matrizes e sistemas lineares

Uma tabela de números que age sobre vetores

Uma matriz é uma tabela retangular de números. Em Python a guardamos como uma lista de linhas, então [[1, 2], [3, 4]] é uma matriz 2×2 e m[linha][coluna] lê uma entrada. Matrizes descrevem muitas linhas do mesmo tipo de cálculo de uma vez. Se cada linha tem os pesos de uma carteira e o vetor tem os retornos do dia dos ativos, o produto da matriz pelo vetor dá o retorno de cada carteira em um passo. Cada entrada do resultado é o produto escalar de uma linha com o vetor, então o comprimento da linha deve ser igual ao do vetor.

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def mat_vec(matrix, vector):
    return [dot(row, vector) for row in matrix]

print(mat_vec([[1, 2], [3, 4]], [5, 6]))   # [17, 39]

# Two portfolios, three assets, one day of returns
weights = [[0.5, 0.3, 0.2],
           [0.2, 0.5, 0.3]]
returns = [0.02, -0.01, 0.03]
print([round(r, 4) for r in mat_vec(weights, returns)])   # [0.013, 0.008]

Transposta e multiplicação de matrizes

A transposta troca linhas e colunas: a entrada (i, j) vira (j, i), então uma matriz 2×3 vira 3×2. Em Python, [list(coluna) for coluna in zip(*m)] faz isso. O produto A @ B de uma matriz n×k e uma k×m é uma matriz n×m cuja entrada (i, j) é o produto escalar da linha i de A com a coluna j de B. As colunas de B são as linhas da sua transposta, o que encurta o código. As dimensões internas devem coincidir (k e k) e, ao contrário dos números, a ordem importa: em geral A @ B não é B @ A.

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def transpose(m):
    return [list(column) for column in zip(*m)]

def mat_mul(a, b):
    columns = transpose(b)
    return [[dot(row, column) for column in columns] for row in a]

print(transpose([[1, 2, 3], [4, 5, 6]]))   # [[1, 4], [2, 5], [3, 6]]

A = [[1, 2], [3, 4]]
B = [[5, 6], [7, 8]]
print(mat_mul(A, B))   # [[19, 22], [43, 50]]
print(mat_mul(B, A))   # [[23, 34], [31, 46]]  (a different answer: order matters)

Resolver um sistema linear

Muitas perguntas têm a forma: achar x tal que A @ x = b. Quanto de cada um de dois instrumentos de hedge anula uma exposição dada? Quais cargas de fatores reproduzem os retornos observados? Isso é um sistema linear. A eliminação de Gauss o resolve transformando as equações num triângulo: para cada coluna escolhe-se como pivô a linha de maior valor absoluto nessa coluna (trocando-a para cima, o que mantém pequenos os erros de arredondamento), subtraem-se múltiplos dela das linhas abaixo para criar zeros e depois se resolve de baixo para cima por substituição regressiva. Se um pivô é (quase) zero, as linhas são dependentes, a matriz é singular e não há resposta única: gere um erro em vez de dividir por zero.

def solve(a, b):
    n = len(a)
    m = [list(map(float, row)) + [float(b[i])] for i, row in enumerate(a)]  # augmented matrix
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(m[r][col]))
        if abs(m[pivot][col]) < 1e-12:
            raise ValueError("the matrix is singular")
        m[col], m[pivot] = m[pivot], m[col]
        for r in range(col + 1, n):
            factor = m[r][col] / m[col][col]
            for c in range(col, n + 1):
                m[r][c] -= factor * m[col][c]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):          # back-substitution
        known = sum(m[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (m[i][n] - known) / m[i][i]
    return x

print(solve([[2, 1], [1, 3]], [3, 5]))                    # [0.8, 1.4]
print(solve([[1, 1, 1], [0, 2, 5], [2, 5, -1]], [6, -4, 27]))  # [5.0, 3.0, -2.0]

try:
    solve([[1, 2], [2, 4]], [1, 2])         # the second row is twice the first
except ValueError as error:
    print(error)                            # the matrix is singular

import math

def dot(a, b):
    return sum(x * y for x, y in zip(a, b))

def mat_vec(matrix, vector):
    return [dot(row, vector) for row in matrix]

# Check a solution by multiplying back: A x should give b again
A = [[2, 1], [1, 3]]
b = [3, 5]
x = [0.8, 1.4]                       # the answer from solve(A, b)
residual = [got - want for got, want in zip(mat_vec(A, x), b)]
print([abs(round(r, 12)) for r in residual])  # [0.0, 0.0]  (tiny rounding errors vanish)
print(all(math.isclose(g, w) for g, w in zip(mat_vec(A, x), b)))   # True

# The identity matrix leaves any vector unchanged
identity = [[1, 0], [0, 1]]
print(mat_vec(identity, [7, 9]))     # [7, 9]

• Uma matriz é uma lista de linhas; mat_vec faz o produto escalar de cada linha com o vetor.
• Transposta com zip(*m); A @ B entrada (i, j) = linha i de A vezes coluna j de B; a ordem importa.
• Resolva A x = b por eliminação com pivô e substituição regressiva; um pivô (quase) zero significa singular.

A seguir: variáveis aleatórias e esperança, a linguagem da incerteza por trás de todo modelo de risco e precificação.

Lab: o que é avaliado

• Aplicar uma matriz a um vetor com mat_vec25 pts
• Trocar linhas e colunas com transpose25 pts
• Multiplicar duas matrizes com mat_mul25 pts
• Resolver A x = b por eliminação com solve25 pts