install.packages("quantmod")

## Installing package into '/home/francesco/R/x86_64-pc-linux-gnu-library/3.5'
## (as 'lib' is unspecified)

library(quantmod)

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: TTR

## Version 0.4-0 included new data defaults. See ?getSymbols.

2.1 Scaricare dati da FRED

Dal sito https://fred.stlouisfed.org/ vogliamo scaricare la serie storica dei tassi T-Bill 3 mesi. Si chiama ‘DTB3’. Scarichiamo DTB3 da FRED in modo che il risultato sia un oggetto di classe zoo (ottimo da usare per serie storiche che contengono date).

getSymbols('DTB3',src='FRED',return.class='zoo')

## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.

## [1] "DTB3"

Per osservare i dati:

View(DTB3)

Per vedere un grafico dei dati:

plot(DTB3)

Per togliere i valori NA (Not Available) dal dataset:

r3=na.omit(DTB3)

Per vedere il grafico della nuova serie non dando nessuna etichetta (label) all’asse x, nessuna etchetta all’asse y e il titolo (main) del grafico:

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')

Per inserire la griglia:

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')
grid()

Calcolo della media e dello scarto quadratico medio (deviazione standard):

mean(r3)

## [1] 4.362548

sd(r3)

## [1] 3.104091

Le statistiche dscrittive dei dati si ottengono con il comando seguente:

summary(r3)

##      Index                 DTB3       
##  Min.   :1954-01-04   Min.   :-0.020  
##  1st Qu.:1970-03-23   1st Qu.: 2.160  
##  Median :1986-07-01   Median : 4.330  
##  Mean   :1986-06-20   Mean   : 4.363  
##  3rd Qu.:2002-09-16   3rd Qu.: 5.910  
##  Max.   :2018-11-29   Max.   :17.140

2.2 Inserimento di linee sui grafici

Per inserire delle rette nel grafico il comando è “abline”. Si vuole inserire una linea orizzontale (h) corrispondente alla media, di colore rosso (col=‘red’) e di tipo (line-type = lty) punteggiato (lty=1, linea continua, lty=2, linea tratteggiata, lty=3, linea punteggiata).

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')
grid()
abline(h=mean(r3),col='red',lty=3)

Si inserisce una linea verticale (v) in corrispondenza delle date della crisi sub-prime (le prime tensioni sul mercato interbancario europeo il 9 agosto 2007 e il fallimento di Lehman Brothers il 15 settembre 2008)

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')
grid()
abline(v=as.Date('2007-08-09'),col='blue',lty=3)
abline(v=as.Date('2008-09-15'),col='blue',lty=3)

Lo stesso risultato si può ottenere usando un comando ‘abline’ unico e inserendo come valore di ‘v’ un insieme delle ordinate su cui costruire le linee verticali. Nei calcoli seguenti diamo alle variabili ‘d1’ e ‘d2’ i valori delle date che ci interessano al fine di semplificare la notazione.

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')
grid()
d1=as.Date('2007-08-09')
d2=as.Date('2008-09-15')
abline(v=c(d1,d2),col='blue',lty=3)

2.3 Inserimento di testo nei grafici

In un grafico si può inserire del testo utilizzando il comando ‘text’. Il primo argomento del comando sono le coordinate del punto dove si vuole che il testo compaia. Poi occorre indicare il testo vero e proprio come argomento di ‘labels’. All’interno delle virgolette del testo, l’opzione ‘’ indica che si deve andare a capo (‘n’ sta per ‘new line’). Altri comandi opzionali possono riguardare il colore del testo (con la stessa sintassi usata nei grafici) oppure la posizione (‘pos’) del testo rispetto alle coordinate date (pos=1 indica la posizione sotto alle coordinate, 2 indica a sinistra, 3 sopra e 4 a destra). Se ‘pos’ non viene indicato si intende che il testo deve comparire al centro delle coordinate. Il comando ‘cex’ indica la dimensione del testo. Il valore di default è cex=1. Un numero minore di 1 indica un testo più piccolo del “normale”, mentre un numero maggiore di 1 indica un testo di grande.

plot(r3,xlab='',ylab='',main='Tasso T-Bill 3 mesi')
grid()
abline(v=d2,col='blue',lty=3)
text(d2,15,labels='Crack \n Lehman \n Brothers',pos=1,cex=0.8)

2.4 Caricare dati mediante xls/csv

Da molti siti internet si possono scaricare dati in formato xls oppure csv. Una volta scaricato il file, si può scegliere l’opzione “Importa Dataset” e poi:

• ‘from text (base)…’ per un file di tipo CSV
• ‘from Excel…’ per un file di tipo xls

In entrambi i casi si aprirà una finestra di dialogo nella quale inserire l’indirizzo del file desiderato. Il file verrà letto in automatico e ci sarà la possibilità di scegliere alcune impostazioni. Per esempio si può scegliere se, all’interno del file, le colonne hanno un nome o meno (cioè se gli elementi della prima riga del database devono essere intesi come nomi delle colonne).

2.5 Il differenziale dei tassi lungo-corto

Sul database FRED si trova la serie storica dei differenziali fra i rendimenti dei titoli di Stato a 10 anni e quelli a 2 anni. La serie storica del differenziale ha il nome ‘T10Y2Y’.

getSymbols('T10Y2Y',src='FRED',return.class='zoo')

## [1] "T10Y2Y"

Grafico:

plot(T10Y2Y); grid()

Ora togliamo i valori NA e ricreiamo il grafico (con un titolo):

delta=na.omit(T10Y2Y)
plot(delta,xlab='',ylab='',main='Differenziale tassi a 10 e a 2 anni')
grid()

Aggiungere una linea orizzontale sullo zero e una linea orizzontale sulla media

plot(delta,xlab='',ylab='',main='Differenziale tassi a 10 e a 2 anni')
grid()
abline(h=0,col='red',lty=3)
abline(h=mean(delta),col='blue',lty=3)

2.6 Lavorare sui dati di un array

Utilizzando la variabile “delta” creata nel paragrafo precedente (con il differenziale dei tassi lungo-corto), vediamo come operare sugli elementi di un array.

Per richiamare un valore specifico si può inserire il numero di posizione di quel valore fra parentesi quadre. Per esempio, il primo valore si ottiene con:

delta[1]

## 1976-06-01 
##       0.68

Se l’elemento è di una serie storica e ha una data ad esso associata, viene visualizzata anche la data (come nell’esempio). Per richiamare un gruppo di elementi, si possono utilizzare i due punti come nell’esempio seguente:

delta[1:3]

## 1976-06-01 1976-06-02 1976-06-03 
##       0.68       0.71       0.70

Per trovare quali elementi soddisfano una certa condizione si può usare il comando “which”. Per esempio, se desidero trovare tutti i valori del differenziale pari a zero do il seguente comando:

which(delta==0)

##  [1]  553 1496 3136 3137 3376 5496 5497 5501 5536 5537 5918 7402 7497 7532
## [15] 7538 7542 7546 7711 7712 7724 7738 7742 7748

N.B. Il comando contiene un doppio uguale perché si tratta di una “condizione”. Il singolo uguale è utilizzato solo per attribuire valori alle variabili

Il risultato dà la posizione di tutti gli elementi che sono zero. Se, invece, desidero visualizzare gli elementi che soddisfano tale condizione, allora devo richiamare, dal vettore “delta” gli elementi che occupano le posizioni mostrate dal comando “which”:

delta[which(delta==0)]

## 1978-08-17 1982-06-03 1988-12-29 1988-12-30 1989-12-13 1998-06-05 
##          0          0          0          0          0          0 
## 1998-06-08 1998-06-12 1998-08-03 1998-08-04 2000-02-10 2006-01-20 
##          0          0          0          0          0          0 
## 2006-06-07 2006-07-27 2006-08-04 2006-08-10 2006-08-16 2007-04-13 
##          0          0          0          0          0          0 
## 2007-04-16 2007-05-02 2007-05-22 2007-05-29 2007-06-06 
##          0          0          0          0          0

In questo modo si trovano anche le date per le quali il differenziale è stato nullo.

2.7 Scaricare dati da Yahoo

Su yahoo finance la serie storica dell’indice S&P500 si chiama ^GSPC. La vogliamo scaricare in formato “zoo” e a partire dalla data del 3 gennaio 1950:

getSymbols('^GSPC',src='yahoo',
           return.class='zoo',
           from='1950-01-03')

## 
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
## 
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).

## [1] "GSPC"

Osservando quanto scaricato, si vede che è formato da 6 colonne, ognuna con un nome

View(GSPC)

Per selezionare una colonna particolare del set GSPC vi sono due modi:

• Primo metodo: indicare il nome della variabile e farlo seguire dal simbolo del dollaro e dal nome della colonna
• Secondo metodo: indicare, tra parentesi quadre, quale colonna si vuole selezionare

Primo metodo

SP=GSPC$GSPC.Close

Secondo metodo (si lascia vuoto lo spazio, tra parentesi quadre, dove bisogna prendere tutta la dimensione). Il seguente comando significa: prendere TUTTE le righe, e prendere la quarta colonna:

SP=GSPC[,4]

Per vedere la nuova variabile:

View(SP)

Il comando “which(condizione)” consente per visualizzare la posizione degli elementi che soddisfano la “condizione”. Il seguente comando visualizza quali elementi di SP sono uguali a 20

which(SP==20)

## [1] 191 194

Il risultato è che gli elementi 191 e 194 sono uguali a 20. Per visualizzarli:

SP[191]

## 1950-10-04 
##         20

SP[194]

## 1950-10-09 
##         20

Per ricevere direttamente i valori di SP che soddisfano la condizione si può scrivere:

SP[which(SP==20)]

## 1950-10-04 1950-10-09 
##         20         20

La condizione può essere anche di < o >, come segue:

which(SP<20)

##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [18]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34
##  [35]  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
##  [52]  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68
##  [69]  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
##  [86]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102
## [103] 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
## [120] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
## [137] 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
## [154] 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
## [171] 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
## [188] 188 189 190 192 195 196 197 198 199 202 203 206 207 208 209 210 211
## [205] 212 213 214 215 216 218 219 220 221 222 223 227 228 229 230 231 232
## [222] 233 234 235 236 237 238 239 240 241 242 243 244 246

Si può indicare una doppia condizione con il simbolo &:

SP[which(SP>19 & SP<20)]

## 1950-06-08 1950-06-09 1950-06-12 1950-06-13 1950-06-22 1950-06-23 
##      19.14      19.26      19.40      19.25      19.16      19.14 
## 1950-09-13 1950-09-14 1950-09-15 1950-09-18 1950-09-19 1950-09-20 
##      19.09      19.18      19.29      19.37      19.31      19.21 
## 1950-09-21 1950-09-22 1950-09-25 1950-09-26 1950-09-27 1950-09-28 
##      19.37      19.44      19.42      19.14      19.41      19.42 
## 1950-09-29 1950-10-02 1950-10-03 1950-10-05 1950-10-10 1950-10-11 
##      19.45      19.69      19.66      19.89      19.78      19.86 
## 1950-10-13 1950-10-16 1950-10-17 1950-10-20 1950-10-23 1950-10-26 
##      19.85      19.71      19.89      19.96      19.96      19.61 
## 1950-10-27 1950-10-30 1950-10-31 1950-11-01 1950-11-02 1950-11-03 
##      19.77      19.61      19.53      19.56      19.73      19.85 
## 1950-11-06 1950-11-08 1950-11-09 1950-11-10 1950-11-14 1950-11-15 
##      19.36      19.56      19.79      19.94      19.86      19.82 
## 1950-11-16 1950-11-17 1950-11-20 1950-11-21 1950-11-28 1950-11-29 
##      19.72      19.86      19.93      19.88      19.56      19.37 
## 1950-11-30 1950-12-01 1950-12-05 1950-12-06 1950-12-07 1950-12-08 
##      19.51      19.66      19.31      19.45      19.40      19.40 
## 1950-12-11 1950-12-12 1950-12-13 1950-12-14 1950-12-15 1950-12-18 
##      19.72      19.68      19.67      19.43      19.33      19.85 
## 1950-12-19 1950-12-20 1950-12-21 1950-12-26 
##      19.96      19.97      19.98      19.92

Per vedere le date di un gruppo di elementi il comando è “time”

time(SP[which(SP>19 & SP<20)])

##  [1] "1950-06-08" "1950-06-09" "1950-06-12" "1950-06-13" "1950-06-22"
##  [6] "1950-06-23" "1950-09-13" "1950-09-14" "1950-09-15" "1950-09-18"
## [11] "1950-09-19" "1950-09-20" "1950-09-21" "1950-09-22" "1950-09-25"
## [16] "1950-09-26" "1950-09-27" "1950-09-28" "1950-09-29" "1950-10-02"
## [21] "1950-10-03" "1950-10-05" "1950-10-10" "1950-10-11" "1950-10-13"
## [26] "1950-10-16" "1950-10-17" "1950-10-20" "1950-10-23" "1950-10-26"
## [31] "1950-10-27" "1950-10-30" "1950-10-31" "1950-11-01" "1950-11-02"
## [36] "1950-11-03" "1950-11-06" "1950-11-08" "1950-11-09" "1950-11-10"
## [41] "1950-11-14" "1950-11-15" "1950-11-16" "1950-11-17" "1950-11-20"
## [46] "1950-11-21" "1950-11-28" "1950-11-29" "1950-11-30" "1950-12-01"
## [51] "1950-12-05" "1950-12-06" "1950-12-07" "1950-12-08" "1950-12-11"
## [56] "1950-12-12" "1950-12-13" "1950-12-14" "1950-12-15" "1950-12-18"
## [61] "1950-12-19" "1950-12-20" "1950-12-21" "1950-12-26"

Se si vuole osservare quale elemento corrisponde alla data del 9 agosto 2007:

which(time(SP)==as.Date('2007-08-09'))

## [1] 14493

Per individuare il prezzo corrispondente a tale data:

SP[which(time(SP)==as.Date('2007-08-09'))]

## 2007-08-09 
##    1453.09

Il grafo dello S&P500 si ottiene con (separando due comandi con “;”, R li esegue entrambi):

plot(SP); grid()

I rendimenti logaritmici sono calcolati nel modo seguente:

dlnSP=diff(log(SP))
plot(dlnSP); grid()

Il rendimento medio annuale è (ricordiamo che i dati sono giornalieri e ci sono 250 giorni lavorativi in un anno):

mean(dlnSP)*250

## [1] 0.07366982

Si rappresenta una linea orizzontale punteggiata sulla media

plot(dlnSP); grid()
abline(h=mean(dlnSP),col='green',lwd=3,lty=3)

Si rappresentano due linee verticali (rosse punteggiate) sulla data dell’inizio della bolla dot-com esulla data dello scoppio della bolla dot-com:

plot(dlnSP); grid()
abline(h=mean(dlnSP),col='green',lwd=3,lty=3)
abline(v=as.Date('1995-01-01'),lty=3,col='red',lwd=3)
abline(v=as.Date('2001-01-01'),lty=3,col='red',lwd=3)

Si rappresentano due linee verticali (blu punteggiate) sulla data dello scoppio della bolla sub-prime e sulla data del fallimento Lehman Brothers:

plot(dlnSP); grid()
abline(h=mean(dlnSP),col='green',lwd=3,lty=3)
abline(v=as.Date('1995-01-01'),lty=3,col='red',lwd=3)
abline(v=as.Date('2001-01-01'),lty=3,col='red',lwd=3)
abline(v=as.Date('2007-08-09'),lty=3,col='blue',lwd=3)
abline(v=as.Date('2008-09-15'),lty=3,col='blue',lwd=3)

2.8 Fondere database con date diverse

Un problema tipico di chi scarica dati da dabase diversi è che le date potrebbero non essere omogenee. Esiste una funzione che permette di fondere tra loro due database basandosi sulle date di riferimento. Creo un finto database di quattro date (con quattro valori estratti da una Normale standard):

date1=as.Date(c('2018-11-01','2018-11-02','2018-11-03','2018-11-04'))
numeri1=rnorm(4,0,1)
esempio1=zoo(numeri1,date1)
esempio1

##  2018-11-01  2018-11-02  2018-11-03  2018-11-04 
## -1.59113849  0.52634734 -0.03244364 -0.75739409

Credo un secondo database (altrettanto finto) con date solo parzialmente uguali:

date2=as.Date(c('2018-11-03','2018-11-04','2018-11-05','2018-11-06'))
numeri2=rnorm(4,0,1)
esempio2=zoo(numeri2,date2)
esempio2

## 2018-11-03 2018-11-04 2018-11-05 2018-11-06 
##  1.5912358  0.1490203 -0.2158472 -0.8162095

Per “fondere” i due database si utilizza il comando “merge”, con i seguenti risultati:

esempio=merge(esempio1,esempio2)
esempio

##               esempio1   esempio2
## 2018-11-01 -1.59113849         NA
## 2018-11-02  0.52634734         NA
## 2018-11-03 -0.03244364  1.5912358
## 2018-11-04 -0.75739409  0.1490203
## 2018-11-05          NA -0.2158472
## 2018-11-06          NA -0.8162095

Osserviamo che viene creato un nuovo database che contiene tutte le date possibili e contiene anche i valori NA in corrispondenza delle date per cui non esistono alcuni valori. Per togliere i valori NA basta utilizzare il corrispondente comando.

esempio=na.omit(esempio)
esempio

##               esempio1  esempio2
## 2018-11-03 -0.03244364 1.5912358
## 2018-11-04 -0.75739409 0.1490203
## attr(,"na.action")
## [1] 5 6 1 2
## attr(,"class")
## [1] "omit"

t=seq(1,length(SP))

reg=lm(log(SP) ~ t)

summary(reg)

## 
## Call:
## lm(formula = log(SP) ~ t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71971 -0.18240  0.00727  0.19581  0.75939 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.100e+00  4.236e-03   731.7   <2e-16 ***
## t           2.747e-04  4.231e-07   649.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2789 on 17340 degrees of freedom
## Multiple R-squared:  0.9605, Adjusted R-squared:  0.9605 
## F-statistic: 4.217e+05 on 1 and 17340 DF,  p-value: < 2.2e-16

a=reg$coefficients[1]
b=reg$coefficients[2]
a;b

## (Intercept) 
##    3.099571

##            t 
## 0.0002747311

plot(SP)

plot(SP)
lines(time(SP),exp(a+b*t),col='red')
grid()

3.1 Regressione su un sottoinsieme

Per selezionare un sottoinsieme della serie di dati si usa il comando ‘window’ che prende come argomenti la serie dei dati, l’eventuale data di inizio della “finestra” e l’eventuale data di fine. Nel nostro caso vogliamo calcolare una regressione dall’inizio della serie fino al primo gennaio 1995:

SP_sub=window(SP,end='1995-01-01')

Si calcola, quindi, la regressione sul sottoinsieme

t_sub=seq(1,length(SP_sub))
reg_sub=lm(log(SP_sub) ~ t_sub)
summary(reg_sub)

## 
## Call:
## lm(formula = log(SP_sub) ~ t_sub)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58953 -0.21342  0.06178  0.21180  0.33956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.243e+00  4.367e-03   742.6   <2e-16 ***
## t_sub       2.383e-04  6.683e-07   356.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2323 on 11317 degrees of freedom
## Multiple R-squared:  0.9183, Adjusted R-squared:  0.9183 
## F-statistic: 1.272e+05 on 1 and 11317 DF,  p-value: < 2.2e-16

a_sub=reg_sub$coefficients[1]
b_sub=reg_sub$coefficients[2]

Adesso possiamo aggiungere al grafico precedente una seconda regressione (per esempio in blu) utilizzando i coefficienti della nuova regressione (le ascisse sono le date dell’indice SP)

plot(SP)
lines(time(SP),exp(a+b*t),col='red')
lines(time(SP),exp(a_sub+b_sub*t),col='blue')
grid()

Una legenda si può inserrire nel grafico attraverso il comando “legend” che prende come argomenti:

• le coordinae di dove si vuol far comparire la legenda; esiste anche la possibilità di inserire una indicazione testuale del seguente tipo: ‘topleft’, ‘topright’, ‘bottomleft’, ‘bottomright’;
• il testo della legenda inserito mediante l’opzione ‘legend’
• tutte le opzioni già viste per i grafici (‘lty’ per il tipo di linea da creare, ‘col’ per il colore)
• ‘bty’ per indicare il tipo di box all’interno del quale si vuole inserire la legenda; per esempio con bty=‘n’ si intende che non si vuole nessuna cornice
• ‘horiz’ è una variabile di tipo T(rue) oppure F(alse) per indicare che si vuole la legenda in formato orizzontale oppure no; il valore di default è ‘T’
• ‘y.intersp’ oppure ‘x.intersp’ (seguiti dall’uguale e da un valore) indicano lo spazio da inserire fra una voce della legenda e l’altra (in verticale ‘y’ oppure in orizzontale ‘x’); il valore 1 indica uno spazio “normale”

plot(SP)
lines(time(SP),exp(a+b*t),col='red')
lines(time(SP),exp(a_sub+b_sub*t),col='blue')
grid()
legend('topleft',col=c(1,'red','blue'),lty=1,
       legend=c('dati','regressione totale','regressione parziale'),
       bty='n')

L’opzione bty=‘n’ indica che la legenda non verrà circondata da una cornice.

X=rnorm(10000)
plot(X)

plot(X,type='l')

pnorm(2)

## [1] 0.9772499

hist(X)

hist(dlnSP)

mu=mean(dlnSP)
sigma=sd(dlnSP)
X=rnorm(length(dlnSP),mean=mu,sd=sigma)

plot(dlnSP); grid()

plot(dlnSP)
lines(time(dlnSP),X,col='green')
grid()

plot(dlnSP)
lines(time(dlnSP),X,col=rgb(0,0,1,alpha=0.3))
grid()

4.1 Simulazione di processi stocastici

Il processo da cui partire è il seguente \[ S(t+dt)-S(t)=S(t) \mu dt + S(t) \sigma dW(t), \] ovvero

\[ S(t+dt)=S(t) + S(t) \mu dt + S(t) \sigma dW(t), \] partendo da un valore dato \(S(t_0)\). I dati necessari sono definiti come segue:

S0=25 # valore iniziale
mu=0.08 # rendimento medio annuale
sigma=0.15 # diffusione annuale
dt=1/250 # intervallo temporale (un giorno, rispetto a un anno)
T=5 # anni

L’equazione può essere risolta partendo da \(S(t_0)\), calcolando il valore successivo e, una volto conosciuto questo, calcolare il valore ancora successivo, attraverso un processo iterativo. Questa procedura è definita “Monte Carlo”. Per iniziare si definisce un vettore “vuoto” di dimensioni pari al numero di giorni da simulare (cioè T/dt):

S=array(NA,dim=c(T/dt,1))

Il vettore “vuoto” può contenere valori di qualsiasi tipo (anche tutti zeri) e, tuttavia, ai fini del controllo dei risultati, può essere utile inserire tutti valori NA (così si verifica facilmente, alla fine del processo, se tutti i valori sono stati effettivamente creati). Il primo valore deve essere quello dato:

S[1]=S0

A questo punto si possono creare i valori successivi con un ciclo for:

for (i in 2:(T/dt)){
  S[i]=S[i-1]+S[i-1]*mu*dt+
    S[i-1]*sigma*rnorm(1,mean=0,sd=sqrt(dt))
}

Si rappresenta graficamente il risultato:

plot(S,type='l'); grid()

4.2 Il caso di N percorsi

Vediamo come creare N percorsi

N=100 # numero simulazioni

L’idea è quella di creare una matrice di simulazioni. Ogni colonna della matrice contiene un percorso simulato, mentre le righe contengono i prezzi di ogni giorno. Lo schema è quello rappresentato di seguito:

Simulazione 1	Simulazione 2	…	Simulazione N
\(S(t_0)\)	\(S(t_0)\)	…	\(S(t_0)\)
\(S(t_1,1)\)	\(S(t_1,2)\)	…	\(S(t_2,N)\)
\(S(t_2,1)\)	\(S(t_2,2)\)	…	\(S(t_2,N)\)
…	…	…	…
\(S(T,1)\)	\(S(T,2)\)	…	\(S(T,N)\)

Su R si crea una matrice “vuota” di valori. Le righe sono i giorni simulati e le colonne contengono i percorsi possibili

S=array(NA,dim=c(T/dt,N))

Si dà lo stesso valore iniziale a tutte le colonne

S[1,]=S0

Si continua con il ciclo for che, adesso, deve creare una riga intera per volta:

for (i in 2:(T/dt)){
  S[i,]=S[i-1,]+S[i-1,]*mu*dt+
    S[i-1,]*sigma*rnorm(N,mean=0,sd=sqrt(dt))
}

Il comando per rappresentare tutte le simulazioni all’interno di una matrice è il seguente

matplot(S,type='l'); grid()

4.3 Media e intervalli di confidenza delle simulazioni

La media delle simulazioni si può calcolare applicando il comando ‘mean’ a tutte le righe di una matrice meidante il comando ‘apply’ che riceve i seguenti argomenti:

• i dati su cui effettuare i calcoli
• la dimensione (righe = 1, colonne = 2) rispetto alla quale effettuare i calcoli
• la funzione da applicare ai dati
• le eventuali opzioni della funzione

Il comando ‘apply’ deve generare una colonna come quella più a destra rappresentata nella seguente tabella:

Simulazione 1	Simulazione 2	…	Simulazione N	Media
\(S(t_0)\)	\(S(t_0)\)	…	\(S(t_0)\)	\(S(t_0)\)
\(S(t_1,1)\)	\(S(t_1,2)\)	…	\(S(t_2,N)\)	\(\frac{1}{N}\sum_{i=1}^{N}S\left(t_{1},i\right)\)
\(S(t_2,1)\)	\(S(t_2,2)\)	…	\(S(t_2,N)\)	\(\frac{1}{N}\sum_{i=1}^{N}S\left(t_{2},i\right)\)
…	…	…	…
\(S(T,1)\)	\(S(T,2)\)	…	\(S(T,N)\)	\(\frac{1}{N}\sum_{i=1}^{N}S\left(T,i\right)\)

I seguenti comandi consentono di:

• rappresentare i prezzi simulati (tutti in grigio)
• aggiungere la curva della media (con spessore 3)

matplot(S,type='l',ylab='',col='lightgray'); grid()
lines(apply(S,1,mean),lwd=3)

Al fine di rappresentare i quantili delle simulazioni, si può utilizzare ancora la funzione ‘quantile’. Se si vogliono osservare il quinto, decimo, novantesimo e novantacinquesimo percentile della prima simulazione, per esempio, si può scrivere il seguente comando:

quantile(S[,1],probs=c(0.05,0.1,0.9,0.95))

##       5%      10%      90%      95% 
## 22.98275 23.65210 27.88320 28.30432

Nel caso delle 100 simulazioni effettuate, dobbiamo applicare la stessa formula ‘quantile’ attraverso il comando ‘apply’. Così, nel nostro caso, possiamo calcolare quinto e novantacinquesimo percentile di tutte le simulazioni come segue.

quant=apply(S,1,quantile,probs=c(0.05,0.95))
dim(quant)

## [1]    2 1250

Vediamo che il risultato è una matrice di due righe, ognuna delle quali contiene la serie dei quantili. La rappresentazione grafica attraverso ‘matplot’, tuttavia, prevede che le serie da disegnare siano contenute nelle colonne (e non nelle righe). Per utilizzare ‘matplot’, dunque, bisogna trasporre la matrice ‘quant’ mediante il comando ‘t()’. I seguenti comandi consentono di:

• rappresentare i prezzi simulati (tutti in grigio)
• aggiungere la matrice trasposta dei quantili (in blu)
• aggiungere la curva della media (con spessore 3)

matplot(S,type='l',col='lightgray',xlab='giorni',ylab='')
matlines(t(quant),type='l',col='blue')
lines(apply(S,1,mean),lwd=3)
grid()

Infine, possiamo inserire una legenda.

matplot(S,type='l',col='lightgray',xlab='giorni',ylab='')
matlines(t(quant),type='l',col='blue')
lines(apply(S,1,mean),lwd=3)
legend('topleft',bty='n',
       legend=c('Simulazioni','5-95 percentili','Media simulazioni'),
       lty=1,lwd=c(1,1,3),col=c('lightgray','blue','black')
)
grid()

Al di fuori del quinto e del novantacinquesimo percentile dovrebbero trovarsi circa 5 percosi rispettivamente.

4.4 Inserimento di aree ombreggiate nei grafici

Il comando ‘polygon’ consente di rappresentare un’area colorata compresa fra due curve. Le coordinate delle curve vanno inserite come i primi due argomenti del comando (prima le ascisse e poi le ordinate). Supponiamo di avere due curve e di voler colorare lo spazio tra di loro. Di seguito, per rappresentare due insieme di valori come se fossero una matrice, li leghiamo in colonna con il comando ‘cbind’ (il corrispondente comando ‘rbind’ unisce due o più insiemi di valori come se fossero righe di una matrice). Nei seguenti comandi viene usato anche il comando ‘rev’ che serve per inveritre i valori di un array (per esempio rev(c(1,2,3)) restituisce i valori 3,2,1).

A=c(5,4,5,4)
B=c(3,2,1,2)
matplot(cbind(A,B),type='l',lty=1,col=1,ylab='')
text(c(1,2,3,4,4,3,2,1),c(A,rev(B)),
     labels=c('A','B','C','D','E','F','G','H'),
     pos=c(1,1,1,1,3,3,3,3))
grid()

Per colorare l’area interna occorre dare al comando ‘polygon’ prima tutte le ascisse dei punti del grafico e poi tutte le ordinate. I punti vanno ordinati in senso orario partendo, per esempio, dal punto A e continuando in ordine alfabetico. Occorre, quindi, dare al comando ‘legend’ tutte le ascisse dei punti da A ad H e poi tutte le ordinate. Poi si può indicare il colore.

matplot(cbind(A,B),type='l',lty=1,col=1,ylab='')
grid()
polygon(c(1,2,3,4,4,3,2,1),c(A,rev(B)),col='blue')

Nella maggior parte dei casi conviene utilizzare un colore trasparente in modo che si possa vedere altre curve o testi sotto all’area colorata. A questo fine bisogna utilizzare il comando ‘rgb(r,g,b,alpha=x)’ dove la sigla indica “Reg-Greeb-Blue” (i colori fondamentali) e gli argomenti sono la percentuale di rosso, di verde e di blu da utilizzare per il colore. Infine ‘alpha’ rappresenta il grado di trasparenza del colore. Se si vuole avere il blu con una trasparenza del 20% il comando è il seguente.

matplot(cbind(A,B),type='l',lty=1,col=1,ylab='')
grid()
polygon(c(1,2,3,4,4,3,2,1),c(A,rev(B)),col=rgb(0,0,1,alpha=0.2))

L’area può essere rappresentata anche senza bordi. Per eliminare i bordi occorre il comando ‘border=NA’.

matplot(cbind(A,B),type='l',lty=1,col=1,ylab='')
grid()
polygon(c(1,2,3,4,4,3,2,1),c(A,rev(B)),col=rgb(0,0,1,alpha=0.2),border=NA)

Nel caso degli intervalli di confidenza, per esempio, si può rappresentare in blu trasparente l’area fra il quinto e il novantacinquesimo percentile. Le ascisse sono date dal numero di elementi contenuti nelle righe dei quantili.

matplot(t(quant),type='l',col='blue',ylab='')
grid()
lines(apply(S,1,mean),lwd=3)
x=dim(quant)[2]
polygon(c(seq(1,x),seq(x,1)),c(quant[2,],rev(quant[1,])),
        col=rgb(0,0,1,alpha=0.2),border=NA)

4.5 Simulazione di processi di Poisson

Per estrarre numeri da una distribuzione di Poisson si utilizza il comando “rpois” che riceve due input: il numero di elementi da simulare e la media del processo (ricordiamo che per un processo di Poisson la media e la varianza coindicono). Creiamo, per esempio, 100 simulazioni, indicando come media (che corrisponde lala “frequenza” dei salti) il valore 0.05. In questo modo ci aspettiamo di osservare circa 5 salti.

rpois(100,0.05)

##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
##  [71] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Per rappresentare graficamente questo processo sipuò utilizzare il comando “plot” che già conosciamo.

P=rpois(100,0.05)
plot(P)

Il processo di Poisson cumulato (ovvero la somma del processo) si ottiene utilizzando la funzione “cumsum” che calcola la somma cumulata dei valori del processo:

cumsum(P)

##   [1] 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3
##  [36] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 5 5 5 5
##  [71] 5 5 5 5 5 6 6 6 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 9 9 9 9 9 9

plot(cumsum(P),type='l'); grid()

Così osserviamo che il processo cumulato di Poisson presenta dei “gradini”.

4.6 Simulare un processo a diffusione e salto

Definiamo un processo di Poisson nella forma differenziale come \(d\Pi(t)\) tale che sia il limite della seguente distribuzione: \[d\Pi\left(t\right)=\begin{cases} 1, & \phi dt\\ 0, & 1-\phi dt \end{cases}\] dove \(\phi\) indica la frequenza dei salti (per ogni unità di tempo \(dt\)). Qui vogliamo simulare un processo nella seguente forma: \[\frac{dS\left(t\right)}{S\left(t\right)}=\left(\mu+\gamma\phi\right)dt+\sigma dW\left(t\right)-\gamma d\Pi\left(t\right)\] dove tutti i parametri sono costanti. Vogliamo rappresentare un fenomeno che: - ha un valore iniziale di 8 euro (\(S(t_0)=8\)) - ha un rendimento medio dell’8% (\(\mu=0.08\)); - ha una volatilità del 20% (\(\sigma=0.2\)); - presenta salti negativi pari al 10%; - presenta frquenza media dei salti di una volta all’anno (cioè \(1\) volta su \(250\) giorni di simulazione).

Se intendiamo simulat \(T=10\) anni per \(N=10\) possibili percorsi, i comandi sono i seguenti:

S0=8
mu=0.08
sigma=0.2
dt=1/250
phi=5
gam=0.1
N=10
T=10
S=array(NA,dim=c(T/dt,N))
S[1,]=S0

for (i in 2:(T/dt)){
  S[i,]=S[i-1,]+
    S[i-1,]*(mu+gam*phi)*dt+
    S[i-1,]*sigma*rnorm(N,0,sqrt(dt))-
    S[i-1,]*gam*rpois(N,phi*dt)
}
matplot(S,type='l'); grid()

Dalle simulazioni i salti potrebbero non sembrare evidenti. Per osservarli meglio si possono rappresentare i rendimenti logaritmici di una qualsiasi delle traiettorie simulate.

plot(diff(log(S[,1])),type='l'); grid()

Da questo grafico si osserva immediatamente che al processo “Normale” sono stati aggiunti dei salti negativi che, in media, sono pari al 10% del valore del titolo.

F=function(x,a){exp(-a*x)-x}

F(10,4)

## [1] -10

F(a=4,x=10)

## [1] -10

F=function(x=0,a=1){exp(-a*x)-x}
F(x=2)

## [1] -1.864665

F(x=2,a=1)

## [1] -1.864665

F(x=c(0,1,2,3))

## [1]  1.0000000 -0.6321206 -1.8646647 -2.9502129

6.1 Grafico di una funzione

Una funzione può essere rappresentata graficamente. Sarà sufficiente creare un array di valori di ascissa e applicare questo array alla funzione.

x=seq(-2,2,0.01)
plot(x,F(x),type='l'); grid()

Ovviamente posso dare alle funzioni un diverso valore del parametro ‘a’ rispetto a quello di default.

plot(x,F(x,a=0.5),type='l'); grid()

6.2 Soluzione di una funzione

Attraverso la funzione ‘uniroot’ si possono ricercare le soluzioni (radici) della funzione ‘F’. Per “radice” si intende il valore dell’argomento della funzione ‘F’ che rende il valore di output uguale a zero. La funzione ‘uniroot’ ricerca il valore di una sola incognita, quindi bisogna dare alla funzione tutti gli altri valori (anche quelli definiti implicitamente). La funzione ‘uniroot’ riceve almeno due input:

• la funzione di cui cercare la radice
• l’intervallo dell’incognita nel quale ricecare la radice
• eventuali altri input, fra cui i valori degli altri input della funzione

uniroot(F,a=1,interval=c(-10,10))

## $root
## [1] 0.5671485
## 
## $f.root
## [1] -8.16764e-06
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

I diversi risultati della funzione si interpretano nel modo seguente:

• $root: si tratta della soluzione ricercata
• $f.root: è il valore della funzione in corrispondenza della radice; ovviamente dovrebbe essere il più vicino possibile a zero
• $iter: numero di iterazioni che sono stati necessarie per raggiungere l’approssimazione desiderata della soluzione; la funzione ‘root’ utilizza un algoritmo che tende a convergere verso la soluzione mediante iterazioni successive; questo metodo è asintotico e converge alla soluzione precisa solo dopo infinite iterazioni, tuttavia la vicinanza alla soluzione è già molto accettabile dopo poche iterazioni
• $init.it: indica di quanto è stato necessario allargare l’intervallo ‘interval’ per trovare la soluzione; se la soluzione si trova dentro l’intervallo indicato all’interno della funzione, allora $init.it restituisce valore NA
• $estim.prec: precisione stimata dell’approssimazione; si tratta della precisione che si ottiene terminando le iterazioni al livello a un numero finito; nel caso si riuscisse ad arrivare alla vera soluzione in un numero finito di iterazioni, $estim.prec assumerebbe valore NA

Può accadere che si dia un intervallo all’interno del quale non esiste la soluzione. In questo caso viene restituito un messaggio di errore. Per verificare si provi a dare il comando:

uniroot(F,a=1,interval=c(1,2))

Si può indicare alla funzione di cercare la soluzione anche all’esterno dell’intervallo dato mediante il seguente comando.

uniroot(F,a=1,interval=c(1,2),extendInt='yes')

## $root
## [1] 0.5671431
## 
## $f.root
## [1] 3.102907e-07
## 
## $iter
## [1] 10
## 
## $init.it
## [1] 6
## 
## $estim.prec
## [1] 6.103516e-05

Osserviamo che, in questo caso, il valore di $init.it non è più NA perché è stato necessario allargare l’intervallo iniziale per trovare la soluzione. Questo ci dà indicazione che l’intervallo era ‘sbagliato’.

6.3 Funzioni con più radici

Quando una funzione ha più radici, il comando ‘uniroot’ deve essere utilizzato con accortezza. Vediamo il caso di una parabola.

F2=function(x){x^2+x-5}

Questa funzione ha due soluzioni reali che possiamo vedere attraverso il seguente grafico.

x=seq(-5,5,0.01)
plot(x,F2(x),type='l');grid();abline(h=0,lty=3)

Se tentiamo di trovare la souzione nell’intervallo (-5,5), la funzione ‘uniroot’ restituisce un messaggio di errore. Si provi a dare il seguente comando:

uniroot(F2,interval=c(-5,5))

Negli estremi dell’intervallo la funzione ha valore positivo e, quindi, ‘uniroot’ suppone che non possa esservi una soluzione. Affinché la funzione valga zero all’interno dell’intervallo, infatti, ‘uniroot’ vuole che essa abbia segno opposto negli estremi. Una soluzione è quella di trovare le due soluzioni separatamente.

uniroot(F2,interval=c(-5,0))

## $root
## [1] -2.791276
## 
## $f.root
## [1] -5.480349e-05
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

uniroot(F2,interval=c(0,5))

## $root
## [1] 1.791281
## 
## $f.root
## [1] -3.269477e-05
## 
## $iter
## [1] 7
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

Un’alternativa consiste nell’utilizzare un pacchetto nuovo.

install.packages('rootSolve')

## Installing package into '/home/francesco/R/x86_64-pc-linux-gnu-library/3.5'
## (as 'lib' is unspecified)

library(rootSolve)

Questo pacchetto contiene la funzione ‘uniroot.all’ che permette di trovare tutte le soluzione all’interno di un intervallo.

uniroot.all(F2,interval=c(-5,5))

## [1] -2.791288  1.791288

Si tratta delle due soluzioni che erano state trovate in precedenza mediante l’utilizzo della funzione ‘uniroot’ per due volte.

6.4 Funzione per simulare processi stocastici

La procedura di simulazione che abbiamo mostrato nei paragrafi precedenti (basata sulla discretizzazione) è dovuta al matematico svizzero Euler. Vogliamo simulare mediante il metodo di Euler la seguente equazione differenziale:

\[dX(t) = \mu \left(t,X(t)\right) dt + \sigma\left(t,X(t)\right) dW(t) \] dando come input le funzioni \(\mu\left(t,X(t)\right)\) e \(\sigma\left(t,X(t)\right)\), che sono la deriva e la diffusione del processo. Qui possiamo scrivere una funzione che abbia come input:

• la deriva del processo da simulare
• la diffusione del processo da simulare
• l’intervallo di tempo dt
• il numero di anni da simulare (T)
• il valore iniziale del tempo (t0) e del processo (x0)
• il numero di processi da simulare (N)

La procedura seguente genera la funzione utilizzando la funzione ‘eval(parse(text=x))’ che consente di valutare come se fosse un’espressione matematica la stringa contenuta nella variabile ‘x’. La fuznione termina richiamando la variabile ‘x’ che è l’output della funzione.

euler=function(deriva,diffusione,dt,T,t0,x0,N){
  mu=function(t,x) eval(parse(text=deriva))
  sigma=function(t,x) eval(parse(text=diffusione))
  x=array(0,dim=c(T/dt,N))
  x[1,]=rep(x0,N)
  t=t0
  for (i in 2:(T/dt)) {
    x[i, ]=x[i-1,]+mu(t,x[i-1,])*dt+sigma(t,x[i-1,])*rnorm(N)*sqrt(dt)
    t=t+dt
    }
  x
}

Adesso possiamo usare la funzione per creare, per esempio, 2 anni di un processo di Wiener. Il processo di Wiener si ottiene ponendo \(\mu\left(t,X(t)\right)=0\) e \(\sigma\left(t,X(t)\right)=1\). Vogliamo iniziare il processo al tempo t0=0 e dalla posizione x0=0 (facendo 100 simulazioni).

W=euler(deriva='0',diffusione='1',dt=1/250,T=2,t0=0,x0=0,N=100)
matplot(W,type='l',xlab='giorni',
        ylab='Processo di Wiener',col='lightgray')
grid()

Volendo rappresentare 10 anni (T=10) di prezzi per un titolo che parte da x0=25 euro e per cui deriva e diffusione sono \(\mu\left(t,X(t)\right)=0.08X(t)\) e \(\sigma\left(t,X(t)\right)=0.15X(t)\), i comandi sono i seguenti.

X=euler(deriva='0.08*x',diffusione='0.15*x',dt=1/250,T=10,t0=0,x0=25,N=100)
matplot(X,type='l',xlab='giorni',
        ylab='Prezzi',col='lightgray')
lines(apply(X,1,mean))
grid()

7.1 Il modello di Vasicek

Uno dei processo stocastici utilizzati per descrivere l’andamento del tasso privo di rischio \(r(t)\) è il seguente (detto di Vasicek dal nome dell’autore che lo ha utilizzato per primo): \[ dr(t)=a(b-r(t))dt+\sigma dW(t)\] Il processo ha le seguenti proprietà:

• il tasso di interesse ha una distribuzione normale (quindi sono amessi tassi di interesse negativi);
• il processo torna sulla media di lungo periodo (rappresentata dal parametro \(b\));
• il processo torna sulla media con una forza data dal parametro \(a\);
• \(\sigma\) misura la volatilità costante del processo.

La stima del parametro \(\sigma\) si ottiene immediatamente dalla formula della varianza del processo: \[\mathbb{V}_{t}\left[dr\left(t\right)\right]=\sigma^{2}dt\] da cui si ottiene la stima: \[\sigma=\frac{\sqrt{\mathbb{V}_{t}\left[dr\left(t\right)\right]}}{\sqrt{dt}}\] Una volta stimato il parametro \(\sigma\), si può discretizzare il processo per stimare gli altri parametri. Il primo passaggio è \[r\left(t+dt\right)=r\left(t\right)+a\left(b-r\left(t\right)\right)dt+\sigma dW\left(t\right)\] e successivamnete si può scrivere \[r_{i}=r_{i-1}+a\left(b-r_{i-1}\right)dt+\varepsilon_{i-1}\] ovvero \[r_{i}=a\cdot b\cdot dt+\left(1-a\cdot dt\right)r_{i-1}+\varepsilon_{i-1}\] Dobbiamo, dunque, stimare la regressione seguente: \[r_{i}=\beta_{0}+\beta_{1}r_{i-1}+\varepsilon_{i-1}\] Quando \(\beta_0\) e \(\beta_1\) sono stati stimati, i valori dei parametri \(a\) e \(b\) sono ottenuti risolvendo il seguente sistema:

\begin{cases} {0}=abdt\ {1}=1-adt \end{cases}

da cui

\[\begin{cases} b=\frac{\beta_{0}}{a\cdot dt}\\ a=\frac{1-\beta_{1}}{dt} \end{cases}\]

7.2 La stima dei parametri sui T-Bill a 3 mesi

Si possono scaricare i dati del tasso di interesse sui T-Bill a 3 mesi dal database FRED (si utilizza il comando as.numeric per trattare i dati come numeri e non come una serie storica).

getSymbols('DTB3',src='FRED')

## [1] "DTB3"

r=as.numeric(na.omit(DTB3)/100)
dt=1/250

La stima di \(\sigma\) si effettua come segue:

sigma=sd(diff(r))/sqrt(dt)
sigma

## [1] 0.01388824

La stima degli altri due parametri si ottiene mediante il comando lm già visto in precedenza. La variabile dipenente è il tasso di interesse dalla seconda osservazione fino all’ultima, mentre la variabile indipendente è lo stesso tasso dalla prima alla penultima osservazione.

regr=lm(r[2:length(r)] ~ r[1:(length(r)-1)])
summary(regr)

## 
## Call:
## lm(formula = r[2:length(r)] ~ r[1:(length(r) - 1)])
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0126513 -0.0002041 -0.0000085  0.0002019  0.0134423 
## 
## Coefficients:
##                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)          1.878e-05  1.190e-05    1.578    0.114    
## r[1:(length(r) - 1)] 9.996e-01  2.222e-04 4499.001   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0008783 on 16217 degrees of freedom
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9992 
## F-statistic: 2.024e+07 on 1 and 16217 DF,  p-value: < 2.2e-16

Vediamo che l’intercetta non appare molto significatica mentre risulta estremamente significativo il coefficiente dell’auto-regressione. Per ricavare i coefficienti e, da questi, i valori dei parametri del modello, servono i seguenti comandi.

beta0=regr$coefficients[[1]]
beta1=regr$coefficients[[2]]
a=(1-beta1)/dt
a

## [1] 0.1041068

b=beta0/(a*dt)
b

## [1] 0.04509253

La doppia parentesi quadra indica che bisogna prendere solo il valore numerico del regressione, senza considerarne l’etichetta

Osserviamo così che il processo tende a tornare su un valore di equilibrio (\(b\)) veramente prossimo alla media del fenomeno:

mean(r)

## [1] 0.04362548

La forza di ritorno sulla media implica un semivita pari a:

log(2)/a

## [1] 6.658038

implicando un processo che impiega un tempo piuttosto lungo per ritornare sul suo valore medio di equilibrio.

7.3 Simulazione e confronto con i valori reali

La simulazione di un processo di Vasicek si può effettuare utilizzando la seguente funzione.

Vasicek=function(r0,a,b,sigma,dt,T,N){
  r=array(NA,dim=c(T/dt,N))
  r[1,]=r0
  for (i in 2:(T/dt)){
    r[i,]=r[i-1,]+a*(b-r[i-1,])*dt+sigma*rnorm(N,0,sqrt(dt))
  }
  r
}

Ora possiamo simulare, per esempio, 10 percorsi con i valori dei parametri che sono stati ottenuti dalla simulazione e confrontarli con i valori veri del tasso di interesse. Ai fini del confronto si possono iniziare le simulazioni proprio dal primo tasso.

Nel codice seguente la funzione floor (che arrotonda al numero intero più vicino per difetto) è utilizzata per avere sempre un numero intero di osservazioni.

r_sym=Vasicek(r0=r[1],a=a,b=b,sigma=sigma,dt=dt,
              T=floor(length(r)*dt),N=10)
matplot(r_sym,type='l',col='lightgray',
        ylim=c(min(r_sym,r),max(r_sym,r)))
lines(r); grid()

Dal grafico osserviamo che le simulazioni non riescono a catturare i momenti di elevata volatilità del tasso e consentono tassi negativi troppo di frequente.

7.4 Il modello di Cox, Ingersoll e Ross

IL modello detto “CIR” ha la forma seguente: \[ dr(t)=a(b-r(t))dt+\sigma\sqrt{r(t)} dW(t)\] dove l’errore non è più omoschedastico (dipende, infatti, dallo stesso tasso di interesse). Al fine di poter effettuare una stima, tuttavia, occorre ricondurre il modello a una forma omoschedastica. In questo caso la trasformazione seguente \[y(r)=2\sqrt{r}\] rende il processo omoschedastico. Applicando il lemma di Ito, infatti, si ottiene \[dy=\left(-\sqrt{r}a+\left(ab-\frac{1}{4}\sigma^{2}\right)\frac{1}{\sqrt{r}}\right)dt+\sigma dW\] ovvero \[dy=\left(-\frac{a}{2}y+\left(2ab-\frac{1}{2}\sigma^{2}\right)\frac{1}{y}\right)dt+\sigma dW\] La stima del parametro \(\sigma\), dunque, si ottiene calcolando la varianza di \(dy\): \[\mathbb{V}_{t}\left[dy\right]=\sigma^{2}dt\] ovvero \[\sigma=\frac{\sqrt{\mathbb{V}_{t}\left[dy\right]}}{\sqrt{dt}}\] Una volta stimato questo parametro, si può effettuare una regressione discretizzando \(dy\): \[y_{i}=y_{i-1}+\left(-\frac{a}{2}y_{i-1}+\left(2ab-\frac{1}{2}\sigma^{2}\right)\frac{1}{y_{i-1}}\right)dt+\varepsilon_{i-1}\] La regressione da effettuare, dunque, è \[y_{i}=\beta_{1}y_{i-1}+\beta_{2}\frac{1}{y_{i-1}}+\varepsilon_{i-1}\] dove

\begin{cases} {1}=1-dt\ {2}=(2ab-^{2})dt \end{cases}

La soluzione di questo sistema è

\[\begin{cases} a=2\frac{1-\beta_{1}}{dt}\\ b=\frac{1}{2a}\left(\frac{\beta_{2}}{dt}+\frac{1}{2}\sigma^{2}\right) \end{cases}\]

Procediamo con i comandi in R.

7.5 La stima dei parametri

Il problema del modello CIR è che non accetta tassi di interesse negativi. Possiamo, dunque, limitare l’analisi ai tassi di interesse pre-crisi sub-prime (cioè prendere i dati fino al 9-8-2007).

r=window(DTB3/100,end=as.Date('2007-08-09'))
r=as.numeric(na.omit(r))
plot(r,type='l'); grid()

Adesso possiamo creare la nuova variabile \(y=2\sqrt{r}\) e calcolarne la volatilità per la stima di \(\sigma\):

y=2*sqrt(r)
sigma=sd(diff(y))/sqrt(dt)
sigma

## [1] 0.05427393

Passiamo ora al modello di regressione lineare per la stima degli altri parametri (\(a\) e \(b\)).

Nella formula di regressione del comando lm dobbiamo inserire anche il termine “-1” per indicare che non vogliamo l’intercetta.

# Variabile dipendente
Y=y[2:length(r)] 
# Prima e seconda variabile indipendente
X1=y[1:(length(r)-1)]
X2=y[1:(length(r)-1)]^(-1)

regr=lm(Y ~ X1 + X2 -1)

summary(regr)

## 
## Call:
## lm(formula = Y ~ X1 + X2 - 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.032091 -0.001240 -0.000025  0.001227  0.034318 
## 
## Coefficients:
##     Estimate Std. Error  t value Pr(>|t|)    
## X1 9.998e-01  1.175e-04 8511.726   <2e-16 ***
## X2 3.887e-05  2.018e-05    1.926   0.0541 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003432 on 13388 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 1.174e+08 on 2 and 13388 DF,  p-value: < 2.2e-16

Adesso possiamo calcolare il valore dei parametri rimanenti.

beta1=regr$coefficients[[1]]
beta2=regr$coefficients[[2]]

a=2*(1-beta1)/dt
a

## [1] 0.09624609

b=(beta2/dt+0.5*sigma^2)/(2*a)
b

## [1] 0.05813963

La semi-vita di questo processo è:

log(2)/a

## [1] 7.201821

ed esso tende a ritornare su un valore di equilibrio che è prossimo alla media del fenomeno:

mean(r)

## [1] 0.05165961

7.6 Simulazione e confronto con i valori reali

La simulazione di un processo CIR si può effettuare utilizzando la seguente funzione.

CIR=function(r0,a,b,sigma,dt,T,N){
  r=array(NA,dim=c(T/dt,N))
  r[1,]=r0
  for (i in 2:(T/dt)){
    r[i,]=r[i-1,]+a*(b-r[i-1,])*dt+sigma*sqrt(r[i-1,])*rnorm(N,0,sqrt(dt))
  }
  r
}

Ora possiamo simulare, per esempio, 10 percorsi con i valori dei parametri che sono stati ottenuti dalla simulazione e confrontarli con i valori veri del tasso di interesse. Ai fini del confronto si possono iniziare le simulazioni proprio dal primo tasso.

r_sym=CIR(r0=r[1],a=a,b=b,sigma=sigma,dt=dt,
              T=floor(length(r)*dt),N=10)
matplot(r_sym,type='l',col='lightgray',
        ylim=c(min(r_sym,r),max(r_sym,r)))
lines(r); grid()

Dal grafico osserviamo che il comportamento delle simulazioni è più simili al percorso storico dei tassi di interesse. Nei momenti di tassi elevati anche la volatilità diviene più elevata e i tassi di interesse raggiungono anche nelle simulazioni i picchi storici.