Appunti di teoria di Statistica II

Lezione 1 - 05/10/2021

Introduzione ai modelli e metodi statistici per dati nel tempo

Quello che andremo a studiare riguarda la previsione, che può essere:

1. Statica

La previsione è statica quando il fenomeno non evolve nel tempo o la previsione non dipende dal passato, ma solo dall'istante temporale in cui interessa fare la previsione. Esempio: modello di regressione dove y è la variabile risposta e ε è l'errore. Con n=3 non si ottengono risultati molto significativi (le bande blu sono molto ampie). Più n è grande, più le bande si restringono. L’obiettivo è prevedere la y in corrispondenza di una x non osservata. α e β sono stimati, ε è l'errore.

2. Dinamica

La previsione è dinamica quando riguarda il futuro e si basa sulla serie storica osservata. Si osserva la y nel tempo, a diversi istanti. Si ha una serie storica e si vuole prevedere cosa succede nell’ultimo istante (ad esempio quando t = 5). Studiando l’andamento dei dati si trova un algoritmo che fa una previsione.

Studieremo i modelli statistici per osservazioni che evolvono nel tempo. Ciò in genere implica che le osservazioni non siano indipendenti come accade nel campionamento della statistica classica ma le osservazioni successive sono (auto)correlate (cioè quanto l’osservazione al tempo t è correlata con le osservazioni precedenti). Per descrivere questo tipo di dati parleremo di modelli per le serie storiche e processi stocastici.

Lo strumento del calcolo delle probabilità che gestisce un insieme di osservazioni nel tempo auto correlate sono i processi stocastici, una famiglia di infinite variabili casuali.

Indice degli argomenti

• Richiami di inferenza statistica
• Normalità e test di normalità
• Metodi di stima: Minimi quadrati e massima verosimiglianza
• Simulazione stocastica: Monte Carlo e Bootstrap
• Modelli di regressione: Regressione multipla e minimi quadrati; Dati di conteggio e GLM; Regressione non lineare e Splines
• Modelli per serie storiche: Autocorrelazione, Stazionarietà, Modelli ARMA e ARIMA, Regressione con errori ARMA, Dati di conteggio
• Modelli parzialmente osservabili e filtro di Kalman

Esempio di serie storica: Random Walk

La y al tempo t1 è pari all'osservazione precedente più uno shock, un’innovazione ε. Questa innovazione è normale e indipendente, con media 0 e varianza sigma quadro e sono tra loro incorrelate.

Nel primo box è rappresentato il cambio giornaliero. Nel secondo box delle variazioni si ha un fenomeno che oscilla intorno a zero, quindi le differenze prime (y_t – y(t-1)) hanno media pari a zero.

Quindi sono distribuite come una normale. Nel secondo grafico la media è zero, la varianza non è molto costante quindi se costruisco l’istogramma dei dati vedo che c’è asimmetria perché la coda sinistra è lunga mentre quella destra è più corta.

È comunque una distribuzione random walk ma non è gaussiana.

Numero di nuovi positivi in Lombardia durante il COVID

È un fenomeno non stazionario. Prima è basso, poi c’è un forte aumento improvviso e poi inizia a calare lentamente, ecc.

AR(2)

AR(2) sta per auto regressivo, di tipo 2. La y al tempo t è uguale a 0.40 per l’osservazione precedente, meno 0.94 per l’osservazione di due istanti fa, più un’innovazione ε. Il processo rimbalza a causa di 0.94 al secondo periodo. È comunque un processo stazionario rispetto alla media. L’osservazione attuale dipende dall’osservazione precedente con coefficiente positivo (minore di uno); dipende anche dall’osservazione ancora prima, quella a due istanti prima e poi c’è un’innovazione. Questa osservazione genera un modello che ha forti oscillazioni intorno allo zero. Descrive un fenomeno oscillatorio.

Qualità dell’aria a Bergamo

D’inverno l’inquinamento è più alto rispetto all’estate. C’è una forte variabilità stagionale. Modello di regressione lineare con errori ARMA e componente stagionale a 12 mesi. Nel tempo il trend è leggermente decrescente.

Temperatura all’interno di una camera di sterilizzazione

Fenomeno non stazionario con oscillazioni guidate dal protocollo di sterilizzazione (tempo in secondi).

Esempio: COVID-19 - Gestione terapia intensiva

Consideriamo il numero di ricoverati in terapia intensiva per COVID-19 in Italia. Obiettivo è mostrare la previsione a un giorno con il modello stimato. Inclusa l’incertezza della previsione. File dati Covid19 dpc Italia.mat, aggiornamento 4 maggio 2020, fonte Dipartimento Protezione Civile.

Inizio slide STAT

Richiami di inferenza e analisi e test di normalità

La distribuzione normale è una distribuzione simmetrica, può avere una certa media o no. I tre grafici rappresentati hanno media zero e varianza 1 (distribuzione normale standard). La media si indica con mu e la varianza con sigma quadro. Quando la x è lontano da mu la densità va a zero con velocità esponenziale elevata al quadrata. La normale va a zero molto in fretta, è simmetrica perché c’è un quadrato. Nel grafico è sempre X, è sbagliato che nel primo c’è Y.

La probabilità che x dista dalla media più o meno di 1 sigma è di circa il 68%. La probabilità che x disti dalla media più o meno di 2 sigma è del 95% circa.

Questi ci aiutano a capire se siamo vicini alla media o se siamo lontani. In base al sigma riusciamo a capire se siamo vicino alla media o più verso le code.

Supponiamo di avere una distribuzione normale perfettamente simmetrica (curva a campana blu). Con la skewness (indice) si misura la simmetria della curva. Il momento terzo (y - mu) è zero perché gli scarti positivi hanno la stessa probabilità degli scarti negativi.

Se invece la curva non è simmetrica poiché ho la coda destra più lunga di quella sinistra (curva verde), allora lo skiù (cioè il momento terzo standardizzato) è maggiore di zero. Se la distribuzione fosse come la curva rossa, lo skewness sarebbe negativo.

sk = skiù = skewness = momento terzo standardizzato (mi dà informazione sul tipo di asimmetria). Se la distribuzione è simmetrica allora sk = 0; se la distribuzione ha la coda destra più lunga allora sk è positivo. Se la distribuzione ha la coda sinistra più lunga allora sk è negativo.

Anche lo skewness stimato è una variabile casuale con distribuzione, in ipotesi di normalità, normale con media zero e una certa varianza.

Altezza delle code

Posso avere due distribuzioni entrambe simmetriche: una Gaussiana con le code che vanno a zero e una con le code più alte (come la T di student). Per capire l’altezza delle code mi aiuta l’indice di curtosi, k. Prendo la variabile standardizzata (cioè con media 0 e varianza 1) e vedo la variabilità che rimane dopo aver tolto la varianza. Standardizzare vuol dire togliere la media e la varianza.

Quello che ha la curtosi più alta è quello che ha le code più alte, a parità di varianza. Si dimostra che se la distribuzione è normale, con qualsiasi media e varianza, allora la normale ha curtosi pari a 3. Quindi se voglio confrontare la mia distribuzione con la normale calcolo la curtosi della mia distribuzione: se k > 3 le code sono più alte, se k < 3 le code sono più basse della normale.

Quella in rosso è la normale, quella verde ha k > 3 perché ha le code più alte (distribuzione più appuntita della normale). Quella blu ha k < 3 perché le sue code sono più basse della normale (distribuzione meno appuntita della normale). L’area totale di ogni curva deve essere 1 e la varianza delle distribuzioni è 1.

Anche la curtosi campionaria è una variabile casuale con distribuzione di probabilità (se il campione proviene dalla normale), per n grande, è approssimativamente normale con media 3 ed una certa varianza che va a zero.

Lezione 2 - 06/10/2021

La t di Student

La t di Student è una distribuzione (di tipo polinomiale) di dati e se ha infiniti gradi di libertà è una normale. Nell’immagine sottostante vediamo che con 2.5 gdl la Curtosi non esiste perché la t di Student non ha tutti i momenti. La Curtosi della t di Student si misura dal 5 grado di libertà in avanti. Al crescere dei gdl la Curtosi diminuisce fino a convergere a 3. Sulla normale avrei una probabilità pari a 0.001, sulla t di student con 300 gradi di libertà (che è circa una normale) è sostanzialmente identica, non cambia molto. Con 300 gdl quindi ho una normale con probabilità 0.001. Quando ho una curtosi di 4, quindi t di student con 10 gdl, la probabilità è triplicata e il tempo di ritorno (cioè il reciproco di questa probabilità, ogni quanti anni in media succede un evento di questo tipo) passa da 945 anni a 323 anni. Se ho 5 gdl, la Curtosi sale a 9 e la probabilità è 5 volte di quella con 300 gdl. Se baso i calcoli sulla Gaussiana ma il mio fenomeno ha una curtosi è 9, sto sbagliando il calcolo del rischio da 1 a 5 (sto sottostimando il rischio di 5 volte). Se invece ho 2.5 gdl con Curtosi infinita poiché l’integrale diverge, la coda è 0.0055. Capire come sono fatte le code è importante se dobbiamo fare l’analisi di rischio. Se dobbiamo invece calcolare l’indice centrale, ossia la media, non è così importante.

Ho un campione casuale semplice, quindi ho n osservazioni con una certa distribuzione che chiamo fn. Questa distribuzione dipende da un parametro θ (può essere un vettore o uno scalare). Quella fn è la distribuzione congiunta di tutti i dati (cioè è il prodotto delle distribuzioni, quindi situazione di indipendenza ed identica distribuzione IID).

Noi avremo una distribuzione congiunta non tenendo conto però dell’indipendenza. In statistica 1 dovevamo stimare la media della popolazione: avevamo la distribuzione della popolazione con media incognita da stimare. Estraggo un campione casuale (x1, … ,xn) da questa distribuzione e volevamo stimare la media tramite la media campionaria:

Qui invece c’è teta che è un generico parametro da stimare tramite θ cappello (funzione dei dati). Per capire l’incertezza di teta cappello è importante capire la distribuzione dello stimatore almeno approssimativamente normale con media = il parametro che voglio stimare, e una certa varianza che chiamo tau quadro fratto n.

L’intervallo di confidenza è la stima più o meno un delta che calcolo prendendo la radice quadrata della varianza t/ e il percentile della normale std. L’obiettivo era stimare la media della popolazione e per farlo si usava la media del campione. Se ora non mi interessa stimare la media, ma un parametro generale θ, ci chiediamo qual è lo stimatore di quel parametro? Un metodo che produce stimatori è il metodo della massima verosimiglianza: se ho una distribuzione di dati, osservo i dati, calcolo la verosimiglianza (cioè la probabilità dei dati) su questi dati tenendoli fissi e facendo variare θ. Infine massimizzo la verosimiglianza rispetto a θ, cioè la stima di θ è quel valore di θ che rende massima la probabilità dei dati che ho osservato. LIK sta per verosimiglianza. La sigla per la stima di massima verosimiglianza è MLE. x1, …, xn sono i dati fissi (FIX) e vado a guardare la verosimiglianza al variare di teta. Quindi la verosimiglianza la calcolo dopo aver osservato i dati. Per i dati ho tanti possibili teta diversi e non so quale è quello giusto. Quello giusto è quello che rende max la verosimiglianza (argmax L(teta)). L’algoritmo di ottimizzazione di L fornisce la deviazione std, nota in via approssimativa. (standard error = deviazione standard)

C’è un teorema (Ralph Lauren) che ci dice che: sotto certe ipotesi, per n tendente all’infinito, θ cappello n è “ottimale” cioè la media di θ cappello n converge al parametro da stimare (quindi teta) ed inoltre la sua varianza è minore o uguale alla varianza di qualsiasi altro stimatore non distorto (quindi ottimale nel senso che ha la varianza più piccola che posso desiderare per n grande).

Verifica d’ipotesi

Ho anche qui un campione, non per forza casuale semplice ma può anche essere auto-correlato. I test parametrici sono test relativi ad un parametro. Ho un’ipotesi H0 che ipotizza un certo valore per il parametro. Mi chiedo: il parametro è uguale ad un certo valore teta 0? Per fare un test di ipotesi serve il livello di significatività α (o rischio di primo tipo). α è la probabilità di rifiutare l’ipotesi nulla ma l’ipotesi nulla è vera. Nel caso in cui lo stimatore θ cappello è gaussiano rifiuto H0 quando i dati e l’ipotesi sono distanti (θ e θ cappello non sono mai uguali). Θ cappello calcolato è il θ cappello dei miei dati, mentre il θ cappello generico è la variabile casuale. Il p-value è la probabilità che lo stimatore disti dalla media più di quanto la stima basata sui miei dati dista dall’ipotesi, sotto H0. Quando ottengo un p-value molto piccolo rifiuto, quando è grande invece accetto.

Supponiamo di prendere la distribuzione della media, sotto H0. Supponiamo che i dati restituiscano un valore di x barra calcolato, cioè il valore che viene dal campione. Il p-value è l’area tratteggiata. Più il p-value è piccolo e più x barra è lontano da mu con zero; più il p-value è grande e più x barra è vicino a mu con zero.

Test non parametrici

I test non parametrici non dipendono da un parametro e ne esistono di diversi tipi. Mi posso chiedere, ad esempio, se i miei dati sono normali. Ho una certa statistica test T che confronta la distribuzione calcolata a partire dai dati con quella ipotizzata. Con la statistica test rifiuto se la T è grande. Alfa è la significatività del test (la qualità dell’errore di primo tipo). Il p-value sarà la probabilità di avere valori più alti di quelli che ho calcolato con i dati.

Test di normalità

Nel test di Jarque e Bera l’ipotesi nulla ci dice che i dati hanno distribuzione normale con una certa media e una certa varianza ignoti. Questo test si basa sulla simmetria e sulla curtosi, sugli indici di skewness e di curtosi. Va a testare che lo skewness sia zero e la curtosi 3. Questo test prende lo skewness al quadrato diviso per la sua varianza (se abbiamo una distribuzione normale lo skewness ha media zero, quindi lo skewness è asintoticamente normale).

Anche la curtosi è asintoticamente normale, con media 3 (quindi k cappello meno tre è una normale, per n tendente all’infinito, con media zero divisa per la sua varianza (n/24)). Quindi prende anche la curtosi al quadrato e la divide per la varianza.

Il test di Bera Jarque prende la statistica JB e la confronta con una chi quadro con due gradi di libertà. Questo test rifiuta per valori alti della statistica JB. Questo test non distingue se la curtosi è alta e bassa perché mette tutto insieme. Se il risultato del test conduce ad una non normalità vuol dire che c’è una certa asimmetria (positiva o negativa) e una certa curtosi (alta o bassa). Se il test rifiuta H0 posso andare a vedere se c’è più skewness o più curtosi.

Questo test è asintotico, cioè la distribuzione che vado a considerare è per n tendente ad infinito quindi conosco la distribuzione. Quindi come si calcola il p-value quando n è piccolo se la distribuzione chi quadro con 2 gdl vale solo per n grande? Sono stati calcolati i p-value tramite simulazione e sono stati tabulati dentro il software (valori approssimati). I parametri da mettere nel comando jbtest sono i dati (Y), alfa e mctol.

Test di adattamento

Per capire se i miei dati sono normali confronto l’istogramma dei miei dati (quello blu) con la distribuzione normale che ha la stessa media e la stessa varianza dei miei dati (curva rossa). Per fare questo confronto raggruppo in classi e vedo le differenze tra le aree calcolate usando l’istogramma e le aree calcolate usando la normale. Prima di tutto devo raggruppare in classi. C1,…,Ck sono k intervalli che mi partizionano l’asse delle ascisse dell’istogramma. Poi calcolo la frequenza assoluta quindi il numero di yt che sono nella j-esima classe, ottenendo le frequenze osservate che confronto con quelle teoriche (probabilità che la variabile casuale y sotto H0 sia proprio nel j-esimo intervallo).

Quindi il test chi quadro di adattamento, a questo punto, va a confrontare le frequenze osservate con quelle teoriche/attese. Se H0 è vera, la media di O è proprio E (quindi le differenze sono piccole). Se H0 è falsa, la media di O è diversa da E e i quadrati li evitano.

Facendo in questo modo la statistica somma dei quadrati diviso le teoriche si approssima con una chi quadrato quindi conosco la distribuzione e posso calcolare il p-value. Per n che tende all’infinito la statistica X quadro è approssimativa.