Este código numpy pode ser vetorizado?

Question

Escrevi a seguinte função para produzir n realizações de um processo CIR para um determinado conjunto de parâmetros:

def cir_simulations(alpha, mu, sigma, delta_t, n, num_sims):

    x = np.reshape(np.array([mu] * num_sims), (-1, 1))

    for i in range(0, n):
        x = np.concatenate((x, np.reshape(x[:, -1], (-1, 1)) + alpha * (
                    mu - np.reshape(x[:, -1], (-1, 1))) * delta_t + sigma * np.sqrt(
            np.reshape(x[:, -1], (-1, 1))) * np.sqrt(delta_t) * np.random.normal(0, 1, size=(num_sims, 1))), axis=1)

    return x

Este código funciona, mas agora estou me perguntando se seria possível remover o loop e vetorizá-lo totalmente. Tive dificuldade em encontrar uma maneira de fazer isso, pois a operação que está sendo executada no loop é recursiva (os valores na próxima coluna da matriz são não linearmente dependentes dos valores das colunas anteriores).

Além disso, esse código poderia ser simplificado? Em particular, sinto que pode haver uma complexidade desnecessária na maneira como estou acessando a última coluna da matriz usando

np.reshape(x[:, -1], (-1, 1))

Answer 1

Como a saída de cada etapa depende da etapa anterior, é impossível eliminar totalmente o loop, a menos que a fórmula possa ser simplificada para permitir um cálculo de várias etapas. No entanto, o código definitivamente ainda pode ser melhorado.

• np.concatenateGeralmente, a chamada é ineficiente, pois requer realocação de memória e cópia de dados. Tanto quanto possível, deve-se pré-alocar toda a memória necessária e atualizá-la iterativamente.

• Alguns cálculos podem ser retirados do loop para melhorar a eficiência. Por exemplo, a np.random.normalchamada pode ser feita apenas uma vez. O mesmo com alpha * delta_t.

• Nomes de variáveis ​​mais significativos podem ser escolhidos.

• É sempre uma boa prática ter docstrings e dicas de tipo em seu código.

Aqui está uma versão de código aprimorado:

def cir_simulations(alpha: float, mu: float, sigma: float, delta_t: float, sim_steps: int, num_sims: int):
    """
    Simulate the CIR process.

    Parameters:
    Input
    ...

    Output
    ...
    """
    output_shape = (num_sims, sim_steps + 1)
    sim_results = np.empty(output_shape)
    sigma_dW = np.random.normal(0, sigma * np.sqrt(delta_t), size=output_shape) 
    alpha_dt = alpha * delta_t
    sim_results[0, :] = r_prev = mu

    for r_t, sigma_dWt in zip(sim_results[1:], sigma_dW):
        r_t[:] = r_prev = (mu - r_prev) * alpha_dt + np.sqrt(r_prev) * sigma_dWt

    return sim_results

Não estou muito familiarizado com o significado de todos os parâmetros da fórmula. Há duas coisas que tenho dúvidas em seus cálculos:

• A taxa \$r\$na fórmula é inicializado com o mesmo parâmetro muusado no cálculo iterativo. É esperado?
• O incremento \$dW_t\$usou uma variação fixa de 1. Esse deveria ser o caso ou deveria ser outro parâmetro de função?