Il Metodo Montecarlo è una tecnica computazionale che si affida a campionamenti casuali ripetuti per ottenere stime numeriche di quantità che sono difficili o impossibili da calcolare analiticamente. L’idea fondamentale è di simulare un processo stocastico (casuale) un gran numero di volte e utilizzare i risultati di queste simulazioni per approssimare il valore desiderato.
È usato in ambiti come:
Finanza (valutazione di opzioni, rischi)
Fisica (simulazioni di particelle)
Statistica (stima di integrali complessi)
Intelligenza Artificiale (reinforcement learning)
Matematicamente:
Formalmente, se vogliamo stimare una quantità [math]\theta[/math], possiamo cercare di trovare una funzione [math]g(X)[/math] di una variabile aleatoria [math]X[/math] tale che [math]\mathbb{E}[g(X)] = \theta[/math]. Generando un campione di realizzazioni indipendenti e identicamente distribuite (i.i.d.) [math]X_1, X_2, …, X_n[/math] dalla distribuzione di [math]X[/math], possiamo stimare [math]\theta[/math] con la media campionaria:
La forza del Metodo Montecarlo risiede nella Legge dei Grandi Numeri. Questa legge afferma che, all’aumentare del numero di esperimenti o simulazioni, la media dei risultati ottenuti si avvicina sempre più al valore atteso teorico. In altre parole, se eseguiamo un gran numero di simulazioni casuali, la media dei risultati sarà una buona approssimazione del vero valore che stiamo cercando di stimare.
La Legge Forte dei Grandi Numeri afferma che se [math]X_1, X_2, …[/math] è una sequenza di variabili aleatorie i.i.d. con media finita [math]\mu = \mathbb{E}[X_1][/math], allora: [math]P\left(\lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} X_i = \mu\right) = 1[/math] Questo significa che la media campionaria converge quasi certamente alla media teorica.
Immagina un quadrato di lato 2r, con un cerchio inscritto di raggio 1. Generando punti casuali nel quadrato, la proporzione di punti che cadono nel cerchio rispetto al totale sarà circa:
L’idea è di utilizzare la probabilità. Se lanciamo un punto a caso all’interno del quadrato, la probabilità che cada all’interno del cerchio inscritto è proporzionale al rapporto tra le loro aree. Ripetendo questo esperimento molte volte, la frazione di punti che cadono nel cerchio ci darà una stima di questo rapporto, da cui possiamo ricavare [math]\pi[/math].
Spiegazione:
L’area del quadrato è [math]r^2 = 4r^2[/math].
L’area del cerchio è [math]\pi r^2[/math].
La probabilità che un punto casuale nel quadrato cada nel cerchio è [math]P(\text{punto nel cerchio}) = \frac{\text{Area cerchio}}{\text{Area quadrato}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}[/math].
Se generiamo [math]N[/math] punti casuali nel quadrato e [math]N_{\text{dentro}}[/math] di essi cadono nel cerchio, allora per la Legge dei Grandi Numeri: [math]\frac{N_{\text{dentro}}}{N} \approx P(\text{punto nel cerchio}) = \frac{\pi}{4}[/math]
Da cui otteniamo la stima di [math]\pi[/math]: [math]\pi \approx 4 \cdot \frac{N_{\text{dentro}}}{N}[/math].
Codice in Python
import random
import math
def stima_pi(num_punti):
dentro_cerchio = 0
for _ in range(num_punti):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
if x**2 + y**2 <= 1:
dentro_cerchio += 1
return 4 * dentro_cerchio / num_punti
print("Stima di π:", stima_pi(1000000))
Stima di π: 3.142056
Il codice genera num_punti coppie di coordinate casuali (x, y) uniformemente distribuite tra -1 e 1 (all’interno del quadrato). Verifica se il punto cade all’interno del cerchio di raggio 1 centrato nell’origine ([math]x^2 + y^2 \le 1[/math]). La proporzione di punti all’interno del cerchio moltiplicata per 4 fornisce una stima di [math]\pi[/math].
Questa è la funzione semicircolare sopra l’asse x, e l’area sotto la curva da 0 a 1 è π/4.
Vogliamo calcolare l’area sotto la curva [math]f(x) = \sqrt{1 – x^2}[/math] tra 0 e 1. Possiamo interpretare questo integrale come il valore atteso di [math]f(X)[/math] dove [math]X[/math] è una variabile aleatoria uniformemente distribuita tra 0 e 1. Moltiplicando questo valore atteso per la lunghezza dell’intervallo (1 – 0 = 1), otteniamo l’integrale.
Sia [math]X \sim U(a, b)[/math] (distribuzione uniforme tra [math]a[/math] e [math]b[/math]). La densità di probabilità di [math]X[/math] è [math]p(x) = \frac{1}{b – a}[/math] per [math]x \in [a, b][/math] e 0 altrimenti.
Il valore atteso di una funzione [math]f(X)[/math] è dato da: [math]\mathbb{E}[f(X)] = \int_{a}^{b} f(x) p(x) , dx = \int_{a}^{b} f(x) \frac{1}{b – a} , dx = \frac{1}{b – a} \int_{a}^{b} f(x) , dx[/math]
Da cui ricaviamo: [math]\int_{a}^{b} f(x) , dx = (b – a) \mathbb{E}[f(X)][/math]
Per stimare l’integrale, generiamo [math]n[/math] punti casuali [math]x_1, x_2, …, x_n[/math] uniformemente distribuiti tra [math]a[/math] e [math]b[/math] e calcoliamo la media dei valori della funzione: [math]\mathbb{E}[f(X)] \approx \frac{1}{n} \sum_{i=1}^{n} f(x_i)[/math]
Quindi, la stima dell’integrale è: [math]\int_{a}^{b} f(x) , dx \approx (b – a) \cdot \frac{1}{n} \sum_{i=1}^{n} f(x_i)[/math].
Codice Python
La funzione integrale_montecarlo implementa esattamente la formula derivata sopra. Genera n numeri casuali uniformemente distribuiti tra a e b, calcola la somma dei valori della funzione f in questi punti, e moltiplica per la lunghezza dell’intervallo (b - a) e divide per il numero di punti n. La funzione integranda [math]f(x) = \sqrt{1 – x^2}[/math] tra 0 e 1 rappresenta un quarto di cerchio di raggio 1, la cui area è [math]\frac{\pi}{4} \approx 0.7854[/math].
import random
import math
def integrale_montecarlo(f, a, b, n):
somma = 0
for _ in range(n):
x = random.uniform(a, b)
somma += f(x)
return (b - a) * somma / n
# Funzione da integrare
f = lambda x: math.sqrt(1 - x**2)
print("Stima integrale:", integrale_montecarlo(f, 0, 1, 1000000))
Valutiamo un’opzione call europea con prezzo strike [math]K[/math], scadenza [math]T[/math], prezzo iniziale del sottostante [math]S_0[/math], tasso privo di rischio [math]r[/math], e volatilità [math]\sigma[/math]. Il prezzo è dato da:
Dove [math]S_T[/math] segue un moto browniano geometrico.
Il prezzo di un’opzione call europea può essere visto come il valore atteso del payoff scontato al tempo presente. Il payoff al momento della scadenza [math]T[/math] è [math]\max(S_T – K, 0)[/math], dove [math]S_T[/math] è il prezzo del sottostante alla scadenza e [math]K[/math] è il prezzo di esercizio. Il prezzo del sottostante si assume che segua un moto browniano geometrico, il che implica che il suo logaritmo segue un moto browniano. Questo permette di simulare i possibili valori di [math]S_T[/math] utilizzando numeri casuali da una distribuzione normale.
Matematicamente:
Il prezzo del sottostante [math]S_T[/math] alla scadenza [math]T[/math] è dato dalla formula: [math]S_T = S_0 \cdot \exp\left(\left(r – \frac{1}{2} \sigma^2\right) T + \sigma \sqrt{T} Z\right)[/math] dove:
[math]S_0[/math] è il prezzo iniziale del sottostante.
[math]r[/math] è il tasso privo di rischio.
[math]\sigma[/math] è la volatilità del sottostante.
[math]Z[/math] è una variabile aleatoria con distribuzione normale standard [math]N(0, 1)[/math].
Il payoff dell’opzione call alla scadenza è [math]\max(S_T – K, 0)[/math].
Il prezzo dell’opzione call al tempo presente è il valore atteso scontato del payoff: [math]C = e^{-rT} \cdot \mathbb{E}[\max(S_T – K, 0)][/math]
Per stimare questo valore atteso con il Metodo Montecarlo, simuliamo un gran numero di volte il valore di [math]S_T[/math] generando numeri casuali [math]Z[/math] da una distribuzione normale standard. Per ogni simulazione, calcoliamo il payoff e poi ne facciamo la media. Infine, scontiamo questa media al tempo presente.
Codice Python
La funzione prezzo_opzione_call simula n_sim possibili valori del prezzo del sottostante alla scadenza [math]S_T[/math] utilizzando un numero casuale Z estratto da una distribuzione normale standard (random.gauss(0, 1)). Calcola il payoff per ogni simulazione, ne fa la media e poi la sconta al tempo presente utilizzando il fattore [math]e^{-rT}[/math].
import random
import math
def prezzo_opzione_call(S0, K, T, r, sigma, n_sim):
payoff_totale = 0
for _ in range(n_sim):
Z = random.gauss(0, 1)
ST = S0 * math.exp((r - 0.5 * sigma**2) * T + sigma * math.sqrt(T) * Z)
payoff = max(ST - K, 0)
payoff_totale += payoff
prezzo_call = math.exp(-r * T) * payoff_totale / n_sim
return prezzo_call
# Parametri
S0 = 100
K = 105
T = 1
r = 0.05
sigma = 0.2
n_sim = 100000
print("Prezzo call europea:", prezzo_opzione_call(S0, K, T, r, sigma, n_sim))
Prezzo call europea: 7.987949940405533
6. Vantaggi del Metodo Montecarlo
Facile da implementare anche per problemi complessi
Si adatta bene ad alta dimensionalità
Ottimo per valutazioni probabilistiche
7. Svantaggi
Richiede molti campioni per ottenere precisione
Lento rispetto a metodi analitici
Dipende fortemente dalla qualità del generatore di numeri casuali
Il problema della Convergenza Lenta
La Legge dei Grandi Numeri: La Garanzia di Convergenza, Non la Velocità
Come accennato, la Legge dei Grandi Numeri ci assicura che, all’aumentare del numero di simulazioni ([math]N[/math]), la nostra stima Montecarlo convergerà al vero valore atteso. Tuttavia, questa legge non ci dice quanto velocemente avverrà questa convergenza. La velocità di convergenza è cruciale per capire quanto sforzo computazionale (quante simulazioni) è necessario per ottenere un risultato con una precisione desiderata.
L’Errore di Stima
In qualsiasi metodo numerico, l’errore è la differenza tra il valore stimato e il valore vero. Nel contesto del Metodo Montecarlo, se stiamo stimando una quantità [math]\theta[/math] con la stima [math]\hat{\theta}_N[/math] basata su [math]N[/math] simulazioni, l’errore è [math]|\hat{\theta}_N – \theta|[/math].
La Varianza dell’Estimatore Montecarlo
La velocità con cui la nostra stima converge è strettamente legata alla varianza dell’estimatore Montecarlo. Consideriamo il caso più semplice in cui stiamo stimando il valore atteso [math]\mu = \mathbb{E}[X][/math] di una variabile aleatoria [math]X[/math] con varianza finita [math]\sigma^2 = \text{Var}(X)[/math]. L’estimatore Montecarlo naturale è la media campionaria delle [math]N[/math] simulazioni indipendenti di [math]X[/math]:
Dove abbiamo utilizzato la proprietà che la varianza di una somma di variabili aleatorie indipendenti è la somma delle loro varianze, e che [math]\text{Var}(X_i) = \sigma^2[/math] per ogni [math]i[/math] poiché le simulazioni sono i.i.d.
Errore Standard
L’errore standard della stima è la deviazione standard dell’estimatore, che è la radice quadrata della sua varianza:
L’errore standard ci fornisce una misura della tipica grandezza dell’errore tra la nostra stima [math]\hat{\mu}_N[/math] e il vero valore [math]\mu[/math]. Il fatto che l’errore standard sia proporzionale a [math]1/\sqrt{N}[/math] significa che la precisione della nostra stima migliora solo con la radice quadrata del numero di simulazioni.
Per ridurre l’errore di un fattore 10, è necessario aumentare il numero di simulazioni di un fattore [math]10^2 = 100[/math].
Per raddoppiare la precisione (ridurre l’errore della metà), è necessario quadruplicare il numero di simulazioni.
Implicazioni Pratiche
Questa lenta velocità di convergenza ha importanti implicazioni pratiche:
Alto Costo Computazionale: Per ottenere stime molto precise con il Metodo Montecarlo, spesso è necessario eseguire un numero enorme di simulazioni. Questo può richiedere una notevole quantità di tempo e risorse computazionali.
Problemi ad Alta Dimensionalità: Sebbene il Metodo Montecarlo si adatti bene a problemi con molte variabili, la lenta convergenza può diventare un collo di bottiglia, specialmente se la varianza della quantità che si sta stimando è elevata.
Confronto con Metodi Deterministici: Molti metodi numerici deterministici (come le regole di quadratura per l’integrazione numerica in una dimensione) possono avere tassi di convergenza molto più rapidi, ad esempio di tipo polinomiale o addirittura esponenziale rispetto al numero di punti utilizzati. Questo è uno dei motivi per cui, per problemi a bassa dimensionalità dove sono applicabili, i metodi deterministici possono essere preferibili al Metodo Montecarlo.
Esempio Intuitivo:
Immagina di voler stimare la media di una popolazione lanciando un dado molte volte. Ogni lancio è una simulazione. La media teorica è 3.5.
Dopo pochi lanci (diciamo 10), la tua media campionaria potrebbe essere abbastanza lontana da 3.5.
Per ottenere una media campionaria molto vicina a 3.5, dovrai lanciare il dado un numero molto grande di volte (centinaia o migliaia).
Questo esempio illustra come, anche in un caso semplice, è necessario un numero elevato di campioni per ottenere una buona approssimazione della media. L’ordine di convergenza [math]1/\sqrt{N}[/math] formalizza questa osservazione.
In sintesi, la convergenza lenta, caratterizzata dall’errore che diminuisce come [math]1/\sqrt{N}[/math], è uno dei principali svantaggi del Metodo Montecarlo e deve essere attentamente considerata nella scelta di questo metodo per la risoluzione di un problema.
8. Varianti e Ottimizzazioni
Monte Carlo Quasi-Casuale (Sobol, Halton): usa sequenze più “uniformi” per migliorare la convergenza.
Importance Sampling: campiona da distribuzioni che enfatizzano le regioni più “importanti”.
Parallelizzazione: perfetto per sistemi distribuiti o GPU.
9. Conclusioni
Il metodo Montecarlo è uno strumento potente e flessibile. Anche se può sembrare semplice, la sua forza risiede proprio nella sua applicabilità generale. È ideale per affrontare situazioni in cui l’incertezza, la complessità o la dimensionalità elevata rendono impraticabili altri approcci.
