Ce code numpy peut-il être vectorisé?

Question

J'ai écrit la fonction suivante pour produire n réalisations d'un processus CIR pour un ensemble donné de paramètres:

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

Ce code fonctionne, mais je me demande maintenant s'il serait possible de supprimer la boucle et de la vectoriser complètement. J'ai eu du mal à trouver un moyen de le faire, car l'opération effectuée dans la boucle est récursive (les valeurs de la colonne suivante de la matrice dépendent de manière non linéaire des valeurs des colonnes précédentes).

De plus, ce code pourrait-il être simplifié? En particulier, j'ai l'impression qu'il peut y avoir une complexité inutile dans la façon dont j'accède à la dernière colonne du tableau en utilisant

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

Answer 1

Puisque la sortie de chaque étape dépend de l'étape précédente, il est impossible d'éliminer entièrement la boucle à moins que la formule ne puisse être simplifiée pour permettre un calcul en plusieurs étapes. Cependant, le code peut encore être amélioré.

• L'appel np.concatenateest généralement inefficace car il nécessite une réallocation de mémoire et une copie de données. Aussi longtemps que possible, il faut préallouer toute la mémoire nécessaire et la mettre à jour de manière itérative.

• Certains calculs peuvent être sortis de la boucle pour améliorer l'efficacité. Par exemple, l' np.random.normalappel ne peut être effectué qu'une seule fois. Même chose avec alpha * delta_t.

• Des noms de variables plus significatifs pourraient être choisis.

• Il est toujours bon d'avoir des docstrings et des indices de type dans votre code.

Voici une version du code amélioré:

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

Je ne connais pas très bien la signification de tous les paramètres de la formule. Il y a deux choses sur lesquelles je doute dans votre calcul:

• Le tarif \$r\$dans la formule est initialisé avec le même paramètre muutilisé dans le calcul itératif. Est-ce attendu?
• L'incrément \$dW_t\$utilisé une variance fixe de 1. Est-ce censé être le cas ou devrait-il s'agir d'un autre paramètre de fonction?