UN PARZIALE PERCORSO DI BASE (1) SULL’ANALISI DI UNA SERIE STORICA REALE, POCO INTUITIVA, COMMENTATO CON IL LINGUAGGIO R E COL MATHEMATICA DI WOLFRAM del dott. Piero Pistoia

Problemi di inquadramento del testo

PRIMA BOZZA DI INDICE A LINKS INTERNI in via di costruzione

Links

1 – PREMESSA sullo stato dell’articolo
2 – IN ANTEPRIMA : la funzione PRDGRAM e l’esercitazione (8 esercizi) sul PERIODOGRAMMA
3 – IL PROLOGO

RIASSUNTO

PARTE Ia

    1. Cenni operativi sui concetti di statistica implicati nell’analisi di una serie storica

    2. Correlogramma ed il Periodogramma

      1. Il Correlogramma ed il Test di Durbin-Watson

      2. Il Periodogramma

    1. Il modello di Regressione Lineare Semplice (RLS)

      1. Prima direzione di ricerca

      2. Seconda direzione di ricerca

        1. Significato dell’analisi dei residui

        2. Stime sulle grandezze della Popolazione

    1. Cenni al significato di media mobile

PARTE 2a

    1. Analisi della serie storica “ Concentrazione Arsenico”

                           Metodo delle “Medie Mobili Centrate” – Modello Additivo

    1. Scopo della ricerca

    2. Analisi preliminare e individuazione di outliers

    3. La serie corretta

    4. Gli Effetti Stagionali e la serie destagionalizzata y1t

    5. Il Ciclo-Trend smussato e la componente casuale

    6. Il modello di regressione lineare semplice e test relativi

      1. Adeguamento del modello di regressione alla popolazione

      2. Il residuo della regressione e l’affidabilità dei tests

 

4 – Cenni al METODO DELLA MEDIA MOBILE
5 – INIZIO AREA FRA PARENTESI

Programmi utili  in R commentati e controllati. Il Correlogramma , la Statistica di Durbin Watson, il Periodogramma (applicato come esercizio a medie trimestrali). Formule trigonometriche delle armoniche costruite dai dati di sfasamento e ampiezza riportati nei risultati.

6 – CENNO A COMANDI DI CALCOLO ED ORGANIZZAZIONE DEI DATI
Filter, matrix e ts di R. Commento sulle prime istruzioni di R (carica dati da file) e processi per automatizzare i ‘conti’
7 – ECCO QUELLO CHE FAREMO CON R: ‘LETTURE’ SUI PROCESSI
8 – INIZIO COPIA SCRIPTS DEL PROGRAMMA CENTRALE
Vari commenti anche difformi e riflessioni anche alternative
9 – PRIMA PARTE IN SINTESI
10 – SECONDA PARTE IN SINTESI
Un altro tentativo sulla caccia ai residui (media mobile 3*3)

11 – L’EPILOGO

EPILOGO

PARTE IIIa

ULTERIORI APPROFONDIMENTI

1 – APPlICHIAMO UNA REGRESSIONE LINEARE MULTIPLA

              1_1 – COME CALCOLARE LA F DI FISHER NELLE RLM ([3] 856-860)

              1_2 – COME CALCOLARE L’ERRORE STANDARD (ES) SUI COEFFICIENTI DI REGRESSIONE NELLA RLM

2 – APPLICHIAMO UNA REGRESSIONE MULTIPLA “PESATA”

3 – AZZARDIAMO UNA PREDIZIONE NEL FUTURO

4 – CONCLUSIONI E SUGGERIMENTI

BIBLIOGRAFIA

12 -APPENDICE1

Il Correlogramma ed il Test di Durbin-Watson – Lettura Correlogramma

13 -APPENDICE2

PROGRAMMI IN BASIC: calcolo Coefficienti di Autocorrelazione, il Test di Durbin-Watson, il Test della                    normale di Lin-Mudholkor, analisi spettrale per il Periodogramma. Calcolo dei coefficienti in una regressione multipla (MLR), calcoli con le matrici, metodo di Cholescki. Calcola il radicando dell’errore Standard delle predizioni con la RLM, calcolo matriciale. Tavole per il Test di Normalità di Lin-Mudholkar e per il Test di Durbin-Watson.

14 -APPENDICE3

Tabelle 1-4 dei risultati sull’analisi della serie storica in studio relative all’articolo “Esempi guidati di statistica applicata” di P. Pistoia

15 -APPENDICE4

Analisi con il linguaggio R della serie storica trimestrale rivisitata e ampliata con periodogrammi risultati e grafici.

16 -APPENDICE5

ARTICOLO PREMESSA: “Il senso comune, l’insegnamento scientifico ed i saperi preposti alle scelte” di P. Pistoia

ARTICOLO COMMENTO: “Analisi di Fourier con commenti su dati reali e simulati con il Mathematica di Wolfram vers. 4.2.” di P. Pistoia

“PROGRAMMI in Mathematica con esercitazioni” di P. Pistoia

Vari esempi analizzati compreso ‘Oscillazione mensile ozono a Montecerboli (Pomarance, Pi), 2007,2011’

L’Esempio 5 si riferisce all’analisi della serie storica concentrazione As detrendizzata.

 

1 – PREMESSA

PREMESSA SULLO STATO DELL’ARTICOLO

Il presente scritto diventa, sempre più articolato ‘nell’andare’, sempre meno lineare, continuando a riempirsi di parentesi, di alternative informatiche, di pause di riflessione, di ritorni e di correzioni (si veda, per es., il caso del periodogramma come function, ormai praticamente risolto, inseribile come modulo all’interno di qualsiasi programma scritto dai lettori) ecc.. Per me è questo il ‘vero’ articolo scientifico col suo ‘travaglio raccontato (trouble)’, denso di stimoli, possibilità nascoste, interferenze casuali… e non lo scritto finale asettico e razionalmente ripulito, che banalizza il percorso. In questa ottica qualcuno ha detto che l’articolo scientifico è un inganno.  Possiamo forse affermare che seguire il ‘processo’  è come un auto-porsi  domande-risposte, attraverso una successione di ipotesi-falsificazioni, una sorta di MAIEUTICA  SOCRATICA che favorirebbe la costruzione del concetto? Il filosofo non insegna nulla ai discepoli, ma piuttosto a scoprire la ‘verità’, che potenzialmente hanno già dentro di loro (per processo co-evolutivo con la Natura), attraverso una successione di argomentazioni su  punti interrogativi. Allora, dal punto di vista educativo-didattico è più importante il percorso o la meta, la storia o l’evento? (meditate, gente, meditate!). Secondo me si apprende molto più e meglio se spingiamo a riflettere sugli errori  rilevati, sulle ipotesi a cammino chiuso, sulle falsificazioni insomma, anche in termini di memoria, che seguire acriticamente un racconto lineare, ‘ripianato’, anche se intrinsecamente coerente. In questa disquisizione aperta si inserisce bene anche l’altro aspetto di un Socrate-docente che, perchè ‘ignorante’,  costruisce insieme al discepolo, senza conoscenze preacquisite (risuonano qui le posizioni di Foerster e Bruner, da richiamare in questo blog).

Per sovrapporre però una ‘lettura’ su video meno discontinua e difficile, che serva come back-ground, una guida all’apprendimento più lineare,  più conforme, meno a ‘frullato di pezzi di concetti’ e quindi forse più facile e più gradevole,  trasferiamo, col titolo ‘IL PROLOGO’, la prima parte dell’articolo originale dello stesso autore (senza l’uso di R, ma di scripts in Qbasic ed Excel), di cui lo scritto in questione voleva essere una ‘lettura rivisitata’ mediata dal linguaggio R e dal Mathematica di Wolfram. Prima delle appendici trasferiamo anche la seconda parte col titolo ‘L’EPILOGO’. L’intenzione è introdurre all’inizio anche un INDICE a link per migliorare l’accesso alle diverse ‘zone mosaico’ dell’articolo. Mi scuso per ‘questo andare’ poco controllato! Se mi rimanesse più energia mentale e ‘tempo di vita’ forse potrei anche rivisitarlo. 

Comunque, un buon apprendistato sarebbe quello di leggere, prima di questo intervento, il primo post dal titolo “Un percorso verso il periodogramma” curato dallo stesso autore. Grazie.

2 – IN ANTEPRIMA

IN ANTEPRIMA

ECCO LA FUNCTION PRDGRAM DEL PERIODOGRAMMA IN R  di Piero Pistoia

FUNZIONE DEL PERIODOGRAMMA in pdf

FUNZIONE DEL PERIODOGRAMMA1 in odt

ATTENZIONE!

Segue una proposta di esercitazione da attivare sulla consolle di R: 1) si incolla la f. PRDGRAM in R e in successione 2) si trasferiscono gli ESERCIZI dell’esercitazione, per es., uno alla volta. Si hanno i dati e grafici in uscita per ogni ESERCIZIO. Ricordarsi, una volta sulla consolle, per prima cosa, sempre azzerare  i dati, che R ha già in memoria, tramite il menù ‘VARIE’ (Rimuovi tutti gli oggetti) e poi introdurre in R, prima di incollare la PRDGRAM, le ‘library’ necessarie (tseries e graphics). 

period_reg_rand0001

Per vedere in odt l’Esercitazione cliccare sotto:

Periodogramma _di_dati_simul trend_random_mod2_3 in odt

0ppure……. continuare a leggere…….


PROPOSTA DI ESERCITAZIONE ANCHE PER FAVORIRE L'ACQUISIZIONE 
INTUITIVA DELLA 'LETTURA' DI UN PERIODOGRAMMA (contenuta nel 
precedente link)

Inizialmente vogliamo simulare ad hoc una serie storica 
'tabellando' n=21 dati da tre funzioni del seno con costante 
additiva 100,con ampiezze rispettivamente 4,3,6 e 'frequenze' 
nell'ordine 2/21, 4/21,5/21 e infine  fasi -pi/2, 0, -1.745, 
con  il comando iniziale di di R: t=c(1:n), usando come base 
per i nostri esempi proprio questa espressione:
 
yt=100+4*sin(2*pi*2*t/n-pi/2)+3*sin(2*pi*4*t/n+0)+
6*sin(2*pi*5*t/n-1.745) #0.

Calcolati i 21 dati yt, attribuendo a t valori da 1 a 21 
nell'espressione precedente, tali dati rappresentano 
proprio lanostra serie storica da sottoporre al 
Periodogramma, una volta precisati i tre valori 
essenziali da passare ad esso (yt,n,m), dove m è il 
numero di armoniche da calcolare; m=n/2-1 se n è 
pari; m=(n+1)/2 se m è dispari. 
Tramite il nostro programma in R calcolammo allora 
i valori di ampiezze e fasi per le prime 10 armoniche 
riscoprendo nei dati le oscillazioni che c'erano.
Per esercizio continuiamo a simulare serie storiche 
modificandol'espressione di base, modificandola anche 
aggiungendo, a scelta, un trend lineare (k*t) e/o 
valori random onde controllare se il Periodogramma 
riesce a"sentire", oltre alle oscillazioni armoniche, 
anche il trend e la componente casuale.
Con l'istruzione '#' elimineremo secondo la necessità 
le linee di programma non utilizzate per lo scopo 
prefissato.
	 
Proviamo, prima, ad applicare il programma su 21 dati 
simulati dalle espressioni di una retta inclinata e da 
una serie random estratta da una distribuzione gaussiana. 
Sceglieremo poi una combinazione di seni interessanti 
più adatta a proseguire l'esercitazione.  
period_reg_rand0002


PERCORSI DA INVESTIGARE
 
par(mfrow=c(1,1))

 #n=21
 #n=240
			
 #t=c(1:n)
 
 # yt=0.5*t # 1
 #si tratta di un ramo di iperbole(?)discendente
 
 #yt=c();yt[1:t]=0
 
 #yt <- rnorm(t,0,1) # 2
 #yt=-4+ 0.5*t + rnorm(t,0,1) # 3
 
#yt=100+4*sin(2*pi*2*t/256-pi/2)+3*sin(4*t/256*2*pi+0)+
6*sin(5*t/256*2*pi-1.745) # 4 
#analisi yt; tenendo come base questa espressione con 
armoniche basse, ro è sulla rampa alta #della 'iperbole' 
e si obnubila il trend.
 
 #yt=100+4*sin(2*pi*2*t/n-pi/2)+3*sin(2*pi*4*t/n+0)+
6*sin(2*pi*5*t/n-1.745) + 0.1*t # 5 
 
#analisi yt_reg
 
 #yt=100+2*sin(2*pi*2*t/n-pi/2)+sin(2*pi*4*t/n+0)+
3*sin(2*pi*5*t/n-1.745) + rnorm(t,0,1)*2 # 6 
 #analisi yt_rnorm: diminuiamo le ampiezze e aumentiamo 
i random
 
 #yt=100+4*sin(2*pi*2*t/n-pi/2)+3*sin(2*pi*4*t/n+0)+
6*sin(2*pi*5*t/n-1.745) + 0.5*t)+(rnorm(t,0,1)-1/2))  # 7 
 #analisi yt_reg_rnorm

 yt <- 6*sin(2*pi*5*t/n)+2*sin(2*pi*30*t/n)+ 
3*sin(2*pi*40*t/n)+0.1*t + rnorm(n,0,1)*2 # 8 

 #questa espressione anche con 'frequenze' alte (30,40) è la 
 #più indicata a dimostrare che il Periodogramma 'scopre' anche trends 
 #e randoms oltre alle oscillazioni sinusoidali.
 
Ora possiamo prevedere che cosa accade se togliamo una 
o due di queste tre,basta far girare il programma nei 
diversi casi. 
 In questo contesto nel prosieguo useremo invece, per 
esercizi, le tecniche di scomposizione di una serie 
storica: proviamo a 'destagionalizzarla' in successione 
con due o tre medie mobili opportune (o magari col 
comando filter di R) per controllare che cosa rimane 
(che cosa accade ai random?). Potevamo anche 
'detrendizzarla prima con una regressione lineare, 
ovvero eliminare i random con una media mobile 3*3 ecc..
period_reg_rand0003

TRACCIA DEI PERCORSI

ESERCIZIO N° 0

n0=256 # può essere cambiato
t=c(1:n0)
yt0=100+4*sin(2*pi*2*t/n0-pi/2)+3*sin(2*pi*4*t/n0+0)+
6*sin(2*pi*5*t/n0-1.745)
yt0 # la serie storica
ts.plot(yt0)
if(n0/2==n0%%2) m0=n0/2-1 else m0=(n0-1)/2
yt0_period=PRDGRAM(yt0,n0,m0)
yt0_period # data in uscita con ampiezza e fase, per il 
controllo
yt0_period$ro # vettore delle ampiezze
ts.plot(yt0_period$ro)

Esercizio N° 1

n01=21
t=c(1:n01)
yt1=0.5*t
yt1 # serie storica
ts.plot(yt1)
if(n01/2==n01%%2) m01=n01/2-1 else m01=(n01-1)/2
yt1_period=PRDGRAM(yt1,n01,m01)
yt1_period #data in uscita comprese ampiezze e fasi
yt1_period$ro #vettore delle ampiezze
ts.plot(yt1_period$ro)

Esercizio N° 2

n2=21 # può essere cambiato
t=c(1:n2)
yt2<- rnorm(t,0,1)
plot(yt2)
yt2 # serie storica
if(n2/2==n2%%2) m2=n2/2-1 else m2=(n2-1)/2
yt2_period=PRDGRAM(yt2,n2,m2)
yt2_period # data in uscita
yt2_period$ro # vettore delle ampiezze
plot(yt2_period$ro)

ESERCIZIO N° 4

n4=256 # può essere cambiato
t=c(1:n4)

yt4=100+4*sin(2*pi*2*t/256-pi/2)+3*sin(2*pi*4*t/256+0)+

6*sin(2*pi*5*t/256-1.745)
yt4 
ts.plot(yt4)
if(n4/2==n4%%2) m4=n4/2-1 else m4=(n4-1)/2
yt4_period=PRDGRAM(yt4,n4,m4)
yt4_period # data in uscita
yt4_period$ro # vettore delle ampiezze
ts.plot(yt4_reg$ro)




ESERCIZIO N° 5

n5=256 # può essere cambiato
t=c(1:n5)

yt5=100+4*sin(2*pi*2*t/256-pi/2)+3*sin(2*pi*2*pi*4*t/256+0)+

6*sin(2*pi*5*t/256-1.745)-0.1*t

plot(yt5,type=”l”)
if(n5/2==n5%%2) m5=n5/2-1 else m5=(n5-1)/2
yt5_reg=PRDGRAM(yt5,n5,m5)
yt5_reg # data in uscita
yt5_reg$ro # vettore delle ampiezze
ts.plot(yt5_reg$ro)
                               ____________________________________________

perio_reg_rand0001ESERCIZIO N° 8
par(mfrow=c(1,2))
n8=100 # può essere cambiato
t=c(1:n8)

yt8=6*sin(5*pi*2*t/n8-pi/2)+2*sin(2*pi*30*t/n8+0)+3*sin(2*pi*40*t/n8-1.745)+rnorm(n8,0,1)*2

ts.plot(yt8)
if(n8/2==n8%%2) m8=n8/2-1 else m8=(n8-1)/2
yt8_reg=PRDGRAM(yt8,n8,m8)
yt8_reg # data in uscita
yt8_reg$ro # vettore delle ampiezze
ts.plot(yt8_reg$ro)

GRAFICO YT8 E PERIODOGRAMMA (Yt8_reg$ro) SENZA IL TREND
period_confronti0001
GRAFICO DI Yt8_reg_rnorm n=240
period_confronti0002

 

GRAFICO Yt8  ANCHE CON IL TREND (serie originale)
 
 period_confronti0004
#RIFLESSIONI
#Se aggiungo il trend 0.1*t a yt8 ottengo il grafico 
precedente. Confrontando il grafico che segue#e quello 
precedente sarebbe interessante approfondire 
intuitivamente perchè col trend le ampiezze
#vengono disturbate tanto più quanto più lentamente 
scende a zero il ramo di 'iperbole'.Sembra #quasi così, 
induttivamente, si possa affermare la regola empirica 
(ipotesi) che armoniche con #frequenze più alte  vengano 
disturbate meno di quelle più basse, che si posizionano 
sul ramo a #pendenza più elevata e con i suoi punti 
più distanti dall'ascissa. Se sommiamo la distanza della 
#base dei picchi dall'asse orizzontale alla cima dei 
picchi l'ampiezza tenderebbe al valore della 
#formula? Se togliamo anche i random da yt8 i tre picchi 
sarebbero poggiati sull'asse orizzontale?#La numerosità 
di yt8 influisce o no sulla velocità con cui si muove 
verso l'asse x la curva del  trend? Cercare di rispondere 
osservando i grafici precedenti.
 period_reg_rand0004

FINE ANTEPRIMA

<A NAME=”punto3″>IL PROLOGO

IL PROLOGO

3 – PROLOGO

COME INTRODUZIONE RIPORTIAMO LA PRIMA PARTE DELLA RICERCA ORIGINALE (SENZA L’USO DI R);  LA SECONDA PARTE VIENE RIPORTATA PRIMA DELLE APPENDICI. 

piero_stat0001

pier_stat0001

pier_statw30001

SE VUOI APPROFONDIRE LE PROBLEMATICHE RELATIVE A FOURIER VEDI L’APPENDIX5

pier_stat0002

pier_stat50001

pier_stat6y0001

pier_stat0005

pier_stat0006

pier_stat90001pier_stat0007

pier_stat0008
pier_stat120001

pier_statz130001

LA COSTRUZIONE SI FA CON L’ANDARE!

 LA FUNCTION DEL PERIODOGRAMMA ora può essere trasferita come modulo in qualsiasi  altro programma scritto da chiunque!  Abbiamo  cercato di correggere  tutti gli scripts dove figurava questa funzione all’interno di questo post.  Vedere di seguito (area definita “fra parentesi”) il funzionamento di  un listato con svariati richiami a questa funzione con proposte di ‘gioco’ con le armoniche su una serie storica reale (serie storica trimestrale) …. Il   listato del periodogramma è lungo e articolato. Nell’analisi di una serie di dati storici con piu’ serie derivate capita spesso di far uso di questo listato per guardare all’interno delle serie. E’ pertanto utile riuscire a scrivere una sola volta questo listato per poi richiamarlo quando serve. Da riorganizzare anche testo e paragrafi. Problemi sorgono anche perché R memorizza all’uscita tutti gli oggetti su cui ha lavorato che tacitamente, pur nascosti, sono ancora disponibili. Questi valori possono interagire sui programmi in via di sviluppo, creando situazioni le più disparate. In generale conviene dal menù ‘varie’ eliminare questi valori prima di far girare o costruire programmi! Si cercherà con calma  di attivare i controlli  anche sugli altri post, dove figura la function PRDGRAM.

ATTENZIONE: I SEGMENTI DELL’ARTICOLO IN GRIGIO CHIARO HANNO UNA BARRA ORIZZONTALE IN FONDO PER MUOVERE LO SCRITTO A DESTRA E SINISTRA, SE LO SCRITTO ESCE DALLO SCHERMO

stat_reg_mlr_blog0001

 FINE PROLOGO

 

               UN PARZIALE PERCORSO DI BASE SULL’ANALISI STATISTICA DI UNA SERIE STORICA REALE POCO INTUITIVA COMMENTATO CON IL LINGUAGGIO R

“Letture” su concetti statistici e su alcuni aspetti della programmazione

Dott. Piero Pistoia

PREMESSA

NB – I GRAFICI OTTENUTI CON IL SUPPORTO DEL PROGRAMMA CORR IN QBASIC (ALLEGATO) E DI EXCEL,  SE RIUSCIAMO A RIDISEGNARLI TUTTI, FACENDO GIRARE GLI SCRIPTS DEL LINGUAGGIO R CHE SEGUONO, QUESTO E’ UN EFFICACE CONTROLLO INTERNO ALLO SCRITTO.

Il file.dati che prenderemo come campione da analizzare riguarda le concentrazioni mensili di arsenico (As) misurate in mg/l nelle acque della Carlina (sorgenti Onore), prov. Siena, nell’intervallo di tempo 1989- 1993 (5 anni, 60 mesi con inizio da gennaio). Dopo interpolazione per i dati mancanti,   un’analisi preliminare (Modello Additivo secondo il Metodo delle Medie Mobili Centrate) porta ad individuare tre residui standardizzati elevati (> 2 in valore assoluto e quindi considerati outliers da eliminare e sostituire con nuova interpolazione,ottenendo così una serie storica corretta, stocastica e discreta; stocastica, nel senso che il futuro è solo parzialmente determinato dai valori del passato e discreta, nel senso che le misure sono fatte in tempi specifici (ogni mese) a uguali intervalli.

Su questa serie (yt=as1) di 60 dati – inserita nel file che si chiama As-Carlina1.csv – e che comunque   verrà esplicitata all’inizio dell’analisi – procediamo “a fare i conti” e a gestirla con R. Questa parte iniziale preliminare verrà trattata successivamente.

Intanto alleghiamo di seguito Il grafico della serie corretta e interpolata (Graf. N.1).

 

priodogramma0001

L’analisi di base di una serie storica procede alla ricerca delle uniformità al suo interno, come TREND, vari tipi di stagionalità periodica (giornaliera, settimanale, mensile, trimestrale ecc.) correlata al carattere dei dati che abbiamo (orari, giornalieri, settimanali,ecc.), cicli con eventuale periodo superiore che esce dal range dei dati (in generale periodo e ampiezza variabili), la componente random, che riassume lo ‘white noise’ ed altro (impulsi erratici). Alleghiamo come informazioni preliminari anche il relativo grafico dell’autocorrelogramma e del periodogramma (GRAF. N. 2, a e b). Si rimanda al loro significato e processo alla Appendice 1 di questo articolo e al Post scritto a nome di P.Pistoia ed altri, facilmente accessibile da questo sito, per es., battendo periodogramma nella finestra ‘Cerca’. Anticipiamo che dal correlogramma (GRAF. N.2 a)  si osservano una stretta convessità intorno al valore 12-13 che supera la fascia dell’errore, una ondulazione dei picchi (forse una oscillazione), un permanere di picchi nella zona positiva (TREND) ed altro e quindi  si evince che i dati della serie al 95% di fiducia, non sono random e dal periodogramma  si nota un picco forse rilevante corrispondente al valore 5  (5 oscillazioni nel range dei dati, cioè 5 oscill. in 5 anni, una oscillazione all’anno, quindi periodo=12 mesi). In dati mensili, una oscillazione periodica di periodo 12 è allora un’ipotesi plausibile.

Scegliamo di procedere, come tentativo, per prima cosa ad eliminare dalla serie storica corretta ( yt o as1) l’oscillazione stagionale prevista dai grafici precedenti. Useremo vari metodi per farlo e confronteremo poi i risultati.

priodogramma0002

 

4 – Cenni al METODO DELLA MEDIA MOBILE

SINTESI SUL METODO DELLA MEDIA MOBILE

Il metodo della media mobile consiste nel sostituire ai valori osservati, valori artificiali corretti, ottenuti effettuando la media di ciascun valore con quelli contigui (per il calcolo vedere, per es.,  [3] pag. 997), ottenendo una nuova serie storica.

Se da una serie storica vogliamo eliminare una oscillazione di un dato periodo, bisogna scegliere, per il calcolo della media, una lunghezza del periodo mobile uguale il più possibile alla lunghezza del periodo dell’oscillazione prevista.

E’ da tener presente che sembra che talora tale metodo abbia il difetto di inserire un ciclo fittizio in una serie storica anche casuale. Abbiamo controllato nel caso della serie trimestrale enucleata da quella in studio (vedere dopo).

Useremo la Media Mobile Centrata di ordine 12 (come suggerito dai grafici preliminari) che di norma elimina l’oscillazione di uguale periodo insieme alle componenti casuali dalla serie originale, trasformando la serie mensile originale (yt o as1,  che inizia con gennaio, APPENDIX3, TABELLA N.1, col.5  ) in una serie storica di dodici termini più corta (la serie Mbt, APPENDIX3, TABELLA N.1, col.6,  che perde i valori dei primi sei mesi e degli ultimi sei, e inizia da luglio). Da porre attenzione che nel processo di scorciamento il primo termine della serie Mbt si riferisce al mese di luglio del primo anno e così via. L’Mbt sottratta da quella originale (as1) ne fornisce una della stessa lunghezza della precedente (48 temini), l’STRD (componente stagionale + random, APPENDIX3, TABELLA N.1, col.7 ), sulla quale operiamo poi per ottenere il Fattore Stagionale costituito da dodici termini, uno per ogni mese (oscillazione in un anno). Per ottenere il Fattore Stagionale corrispondente ad un mese, si considerano tutti i valori della serie STRD (più corta di 12 termini) corrispondenti a quel mese e se ne fa la media. Quando faremo girare il programma scritto con R e vedremo i 48 valori della serie STRD, potremo controllare che, per es., i 4 valori del mese di gennaio (il settimo, il diciannovesimo, il trentunesimo, il quarantaduesimo) sono -0.0030, -0.0046, 0.0033, 0.0126 e facendo la media otterremo il 7° elemento del Fattore Stagionale, 0.0022, cioè il primo elemento di ESAs (APPENDIX3, TABELLA N.2, col.1), EFFETTO STAGIONALE,  la cui oscillazione è visibile nel GRAF. N.3 a.

Così per il mese di gennaio si fa la media dei 4 valori di gennaio contenuti nella serie STRD, ottenendo il primo valore dell’Effetto e così via. Con questi processi di media verranno eliminate anche le componenti casuali, se ci sono rimaste, dalla serie STRD che diviene così ST (stagionalità). Ripetendo 5 volte la ST copriamo i 5 anni, ottenendo l’Effetto Stagionale. E’ necessario però prima riorganizzare i 12 termini del Fattore Stagionale, spostando i primi sei termini, alla fine degli ultimi sei in maniera da avere i 12 valori allineati da gennaio a dicembre. Per il controllo di questa oscillazione applichiamoci, per es., il programma CORR scritto in Qbasic dall’autore (nota 2) o in linguaggio R (vedere sotto PARENTESI) e focalizziamo l’attenzione sul periodogramma dell’ultima serie ottenuta per osservare la frequenza di questa oscillazione (GRAF. N.3 a,b dell’Effetto Stagionale, ottenuto invece per mezzo di Excel): chiaramente significativa appare la frequenza 5.  Troveremo lo stesso periodogramma anche con R.  Con R useremo la funzione acf (file, main=”Titolo”), per ritrovare i correlogrammi costruiti con CORR ed excel; per il periodogramma si rimanda anche alla relativa routine qui riproposta, rivisitata e funzionante.

————————————————-

5 – INIZIO AREA FRA PARENTESI

5-AREA FRA PARENTESI

APERTA PARENTESI

Alcuni programmi in R utili nello studio delle serie storiche

Da notare (fra parentesi) il programmino riportato qui sotto, scritto in linguaggio R dal sottoscritto, con i suoi risultati, che calcola egregiamente (almeno sembra) i coefficienti di auto-correlazione di una serie storica di prova:

y=c((1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)). Comunque, nell’andare, lo vedremo in azione per i tanti confronti e prove! Si aggiungono di seguito anche scripts in R per il calcolo di DW (test di Durbin Watson), metodo più efficace nell’analisi dei correlogrammi, sempre del sottoscritto.

ATTENZIONE!  GLI SCRIPTS DEI PERIODOGRAMMI COME SUBROUTINES (functions) SONO IN VIA DI CORREZIONE

RIPORTIAMO SUBITO ANCHE IL PROGRAMMA PIU’ COMPLESSO PER COSTRUIRE IL PERIODOGRAMMA DI UNA SERIE STORICA con i  relativi risultati per il controllo . Un controllo quantitativo più puntuale è stato condotto col MATHEMATICA 4.2 di Wolfram nella APPENDIX4 (Piero Pistoia)

Queste routines  messe sotto forma di Functions serviranno per costruire correlogrammi, tests di DW e periodogrammi ognivolta che servono.

library(tseries)

# PROGRAMMINO ‘CORRELOGRAMMA’

# Un piccolo strumento per allenare anche l’intuito

#dott. Piero Pistoia

result=c() # result=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
result1=c() # result1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
#y=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
y=c(1:20) 
# Il lettore può a piacere aggiungere altre funzioni (anche numeri casuali), tentare di indovinare # con ipotesi e poi controllare, per acquisire intuizione sul Correlogramma e sui suoi limiti.

#Controllare se le definizioni dei vettori con elementi NA sono necessari! Sembra di no!
#y=c(1,2,3,4,5)
 N=length(y)
 m=10
 yM=mean(y)

 for(h in 1:m){
for (t in 1:N-h){
 result[t]=(y[t]-yM)*(y[t+h]-yM)
 }
result1[h]=sum(result)
} # OK
result1
result2=c()
#result2=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
#for(h in 1:m){

 for(t in 1:N){
 result2[t]=(y[t]-yM)^2
 }
result3=sum(result2)

# Calcolo il coeff. di correl. di lag 1

rh=result1/result3

t=seq(1:10)

Prh=plot(t,rh)

RISULTATI DELLA PROVA (nessun errore rilevato dalla consolle di R nella prima prova!)

> load(“C:\\Users\\Asus\\Documents\\.RData”)
> library(tseries)

‘tseries’ version: 0.10-32

‘tseries’ is a package for time series analysis and computational
finance.

See ‘library(help=”tseries”)’ for details.

> result=c(); result=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
> result1=c(); result1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
> y=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
>
> #y=c(1,2,3,4,5)
> N=length(y)
> m=10
> yM=mean(y)
>
> for(h in 1:m){
+ for (t in 1:N-h){
+ result[t]=(y[t]-yM)*(y[t+h]-yM)
+ }
+ result1[h]=sum(result)
+ } # OK
Ci sono 45 avvisi (usare warnings() per leggerli)
> result1
[1] 565.25 385.75 233.75 107.25 4.25 -77.25 -139.25 -183.75 -212.75
[10] -228.25
>
> result2=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
> #for(h in 1:m){
>
> for(t in 1:N){
+ result2[t]=(y[t]-yM)^2
+ }
> result3=sum(result2)
>
> # Calcolo il coeff. di correl. di lag 1
>
> rh=result1/result3
>
> t=seq(1:10)
>
> Prh=plot(t,rh)

Risultato da confrontare con acf(y)

SE SCRIVIAMO coeffcorr=acf(y), R DARA’ ANCHE IL VETTORE DATI IN coeffcorr

La formula usata è quella senza la moltiplicazione per N/(N-1)

LA STATISTICA DI DURBIN WATSON

library(tseries) 
y=c(1,2,3,4,5,6,7,8,9,10) 
n=length(y) 
#result=c(NA,NA,NA,NA,NA,NA,NA,NA,NA)
 result=c()
 result1=c()
for(t in 2:n){
 result[t]=(y[t]-y[t-1])^2
}
result=result[2:n]
a=sum(result)

for(t in 1:n)
result1[t]=y[t]
b=sum(y)
dw=a/b
dw

#Nella tabella, k'=n° regressori non contando la costante, a=n° osservazioni (y) e dw, sono le tre informazioni per fare il 
test con la tabella.
#Per k'=1 e a=20  l'intervallo dl-du=1.201-1.411, per cui 0.2 < dl:  presenza di correlazione,
#si respinge l'ipotesi nulla (ipot. nulla = i dati non sono 
correlati!), come era intuitivamente già nelle cose.
Da notare che normalmente il test si applica ai residui per 
testare la loro indipendenza.
RISULTATI DELLA PROVA (nessun errore sulla consolle di R) 
> library(tseries) 
> y=c(1,2,3,4,5,6,7,8,9,10) 
> n=length(y) > 
#result=c(NA,NA,NA,NA,NA,NA,NA,NA,NA) 
> result=c() > result1=c() 
> for(t in 2:n){ + result[t]=(y[t]-y[t-1])^2 + } 
> result=result[2:n] 
> a=sum(result) 
> > for(t in 1:n) 
+ result1[t]=y[t] 
> b=sum(y) 
> dw=a/b 
> dw [1]
 0.1636364
 >

#TENTIAMO SCRIPTS del PERIODOGRAMMA IN FORMA DI FUNCTION del dott. Piero Pistoia

# PROVE_TEST SUL PERIODOGRAMMA E CONTROLLO COL MATHEMATICA 4.2 
# Oscillazioni su medietrim e costruzione delle formule 
trigonometriche
# Eliminazione delle varie armoniche

par(ask=T)
par(mfrow=c(1,3))
#medietrim sono i 20 valori trimestrali relativi ai 60 dati mensili delle concentrazioni arsenico 
#della Carlina per 5 anni, in studio.
#Vedere il Post a nome di Pf.Bianchi_P.Pistoia  "Un percorso verso il periodogramma" 
#in questo blog o rivisitato ed esteso in APPENDIX4.

yt=c(0.04233333,0.06100000,0.04500000,0.0556666,0.05400000,
0.06500000,0.07066667,0.04633333,0.05833333,0.06533333,
0.08516667,0.06866667,0.07650000,0.0761666,0.07300000,
0.06700000,0.07966667,0.07333333,
0.07866667,0.06266667)

#ALTRA PROVA IN COSTRUZIONE
#yt= qui si introduce il vettore detrend_trim, cioè i 20 valori di yt detrendizzato, 
#su cui faremo agire la function del periodogramma. Vedere  
APPENDIX4
# detrend_trim=c(-0.0094714286, 0.0077825815, -0.0096300752, 
#-0.0003760652, -0.0034553885, 
# 0.0061319549, 0.0103859649, -0.0153600251, -0.0047726817, 
0.0008146617, 
# 0.0192353383, 0.0013226817, 0.0077433584, 0.0059973684, 
0.0014180451, 
#-0.0059946115, 0.0052593985, -0.0024865915, 0.0014340852, 
-0.0159785714)
 
n=length(yt)
yt=as.vector(yt)
nx=n
yx=yt 
medietrim=yt




#m =(n-1)/2 # perchè n dispari
#m =(n/2-1) # perchè n pari

if (nx/2%%2==2) mx=nx/2-1 else mx=(nx-1)/2 #controllo 
automatico di n (pari o dispari?)
#Controllare se ho invertito le due opzioni!

nx
mx
t=c(1:length(medietrim))
PRDGRAM<- function(y1,n1,m1) {

# VALORI DEL PARAMETRO ak
a0=c(); k=0; a0=0;
for(t in 1:n1){a0=a0+y1[t]*cos(2*pi*t*k/n1)}
a0

a0=a0*2/n1;a0=a0/2

a0

a=c();a[1:m1]=0;
for(k in 1:m1) {
for(t in 1:n1){
a[k]=a[k]+y1[t]*cos(2*pi*t*k/n1)}}
a=2*a/n1

# vALORI DEL PARAMETRO bk

b=c();b[1:m1]=0;b0=0;k=0
for(k in 1:m1) {
for(t in 1:n1){
b[k]=b[k]+y1[t]*sin(2*pi*t*k/n1)}}

a <- as.vector(a)

for(i in 1:m1){
if (abs(a[i]) < 1e-10) a[i]=0 else a[i]=a[i]}
a

for(i in 1:m1){
if (abs(b[i]) < 1e-10) b[i]=0 else b[i]=b[i]}
b=2*b/n1
b
# AMPIEZZE
#ro[1:m1]=0
ro <- sqrt(a^2 +b^2)

for(i in 1:m1){
if (abs(ro[i]) < 1e-10) ro[i]=0 else ro[i]=ro[i]}

# CALCOLO DELLA FASE DI OGNI ARMONICA
# RIPORTANDO IL VALORE AL QUADRANTE GIUSTO
f2=c()
f2[1:m1]=0
for(i in 1:m1){
f2[i] <- abs(a[i]/b[i])
f2[i] <- atan(f2[i])*180/pi}
f2 =as.vector(f2)
f2
#f2[1:m1]=0 un f2[1:m1] di troppo!
phi <- c()
for(i in 1:m1){
# f2 <- abs(a[i]/b[i]);
# f2 <- atan(f2)*180/pi;
if(b[i]>0 & a[i]>0) phi[i] = f2[i];
if(b[i]<0 & a[i]>0) phi[i] = 180-f2[i];
if(b[i]<0 & a[i]<0) phi[i] = 180+f2[i];
if(b[i]>0 & a[i]<0) phi[i] = 360-f2[i];
if(b[i]==0 & a[i]==0) phi[i] = 0;
if((b[i]<0 & b[i]>0) | a[i]==0) phi[i]=0; 
if(b[i]==0 & a[i]>0) phi[i]=90;
if(b[i]==0 & a[i]<0) phi[i]=360-90
}

# PHI FASE ARMONICHE

phi=as.vector(phi)
phi
param_a <-a
param_b <-b
ampiezza <- ro
fase <- phi

a;b;ro;phi
# Qui, al termine della function si pone il valore di un'unica 
# variabile che esce o, se escono più variabili, si usa  
# un data.frame: data=data.frame(x1,x2,...).
# Ogni chiamata alla function permette di includere l'unica 
# variabile o i data nel nome della chiamata:
# es. periodxx=nome.function(x1,x2,...)

data <-data.frame(a,b,ro, phi) 
data
# questa matrice esce dalla function e viene 'raccolta' nella variabile periodxx

}

#FINE SUBROUTINE ANALISI FOURIER

period=PRDGRAM(medietrim,nx,mx)
period 
plot(period$ro,type="l",main="PERIODG.medietrim",
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")
# 1° grafico in A1
# medietrim (vedere ro del  period. di medietrim) presenta 
# le armoniche rilev. n.3 e n.5 (GRAF.A1)

# for(i in 1:10000000) i=i
#data <-data.frame(param_a,param_b,ampiezza, fase)
#data
# Con il numero delle armoniche considerate rilevanti, 
le relative ampiezze e fasi possiamo
# costruire le loro espressioni trigonometriche.

w1=c(1:length(medietrim))
y_osc=0.0058*sin(2*pi*5*t/20+3.9) # questa oscillazione 
dovrebbe avere  
# un'armonica 5 (GRAF.A3)
so=medietrim-y_osc # so nel grafico dell'ampiezza (GRAF.B2). 
# Questa sottrazione eliminerà l'armonica 5
#  da ro di medietrim (GRAF.B2)

so
#PER UN'ALTRA PROVA

# Se consideriamo l'altra espressione y_osc1=0.0066*sin(2*pi*3*t/20+2.92), che ha un picco 
#all'armonica 3, invece di y_osc, e la sottraiamo da medietrim che ha pure un picco  
#all'armonica 3 (GRAF.A1), come diverrà il grafico? (vedere 
GRAF.B3)

#Se detrendiziamo medietrim (detrend_trim) e applichiamo il 
period. 
#potremo controllare le sue armoniche rilevanti e esprimere in forma analitica 
#(in formula trigonometrica) la loro rilevanza (y_oscxx). 
APPENDIX4 

#detrend_trim=c(-0.0094714286, 0.0077825815, -0.0096300752, 
#-0.0003760652, -0.0034553885, 
#0.0061319549, 0.0103859649, -0.0153600251, -0.0047726817, 
0.0008146617, 
#0.0192353383, 0.0013226817, 0.0077433584, 0.0059973684 
0.0014180451, 
#-0.0059946115, 0.0052593985, -0.0024865915, 0.0014340852, 
-0.0159785714) #ripreso dall'APPENDIX4
 
FINE ALTRA PROVA
ny=length(y_osc) 
n=length(so) 

if (n/2== n%%2) m=n/2-1 else m=(n-1)/2 
period1=PRDGRAM(so,n,m) 
period1 
period1$ro 
#plot(period1$ro,type="l",main="PERIODG.senza osc.5", 
#xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")
y_osc1=0.0066*sin(2*pi*3*t/20+2.92)# armonica 3; FIG.A2 
nz=length(y_osc1)
if (nz/2== nz%%2) mz=nz/2-1 else mz=(nz-1)/2
period6=c() 
period6=PRDGRAM(y_osc1,nz,mz) 
period6 
plot(period6$ro,type="l",main="PERIODG.y_osc1",
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")# 2° grafico in A2
if (ny/2== ny%%2) my=ny/2-1 else my=(ny-1)/2 
period2=PRDGRAM(y_osc,ny,my)  
period2  
period2$ro  
plot(period2$ro,type="l",main="PERIODG.y_osc", 
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza") # 3° grafico in A3
 
period3=c() 
period3=period 
plot(period3$ro,type="l",main="PERIOD.medietrim", 
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")
# 4° grafico in B1 
# medietrim (vedere ro del period. di medietrim)

 

so1=medietrim-y_osc1 
#period4=c() 
#period4=period1 
#plot(period4$ro,type="l",main="PERIODG.senza osc.3", 
#xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")

nz=length(y_osc1) 
if (nz/2%%2==2) mz=nz/2-1 else mz=(nz-1)/2 #controllo automatico di n (pari o dispari?) 
period5=c() 
period5=PRDGRAM(so1,nz,mz) 
period5 
plot(period5$ro,type="l",main="PERIODG.senza osc.3",  
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza") # 5° grafico in B3 
#par=(mfrow=c(1,1)) 
#period6=c() 
period6=PRDGRAM(y_osc1,nz,mz) 
#period6 
#plot(period6$ro,type="l",main="PERIODG.y_osc1",  
#xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")# 
plot(period1$ro,type="l",main="PERIODG.senza osc.5", 
xlab="Armoniche = N° oscillazioni in n dati", ylab="ampiezza")
#6° grafico in B2
#RISULTATI OK
cliccare qui sotto per vedere i risultati degli scripts in odt che verranno costruiti facendo girare il programma precedente.
PERIOD_PROVE_TEST  OK                                                                     

Si aggiungono qui i relativi tre grafici FIG.A, FIG.B, FIG.C costruiti dal programma precedente, e la successiva  FIG.D, che illustra, alla rinfusa, l’appunto relativo alla formulazione delle due armoniche costruite su ampiezze e fasi dei risultati.

FIG.A0001

FIG.A0002

FIG.A0004

FIG.D

FIG.D0001

DA QUI IN POI QUALCOSA ANCORA DA CONTROLLARE

PER VEDERE LA PRIMA VERSIONE DEL PRECEDENTE PROGRAMMA IN PDF 
CLICCARE SOTTO: 
funzion_period_OK_3_richiami_result. 
LA NUOVA VERSIONE DEL PRECEDENTE PROGRAMMA CON IN USCITA 12 GRAFICI SI TROVA CLICCANDO SU: 
FUNCTION_PERIOD_OK_3_RICHIAMI_RESULT 
Una volta compreso come richiamare e come gestire i risultati 
della function del periodogramma, 
ora siamo in grado di continuare di volta in volta la 
correzione. 
#In ogni caso gli scripts dei programmi presentati in R possono essere trasferiti, anche 
#un pezzo alla volta, direttamente sulla console di R con Copia-Incolla: il programma inizierà 
#nell'immediato a girare costruendo risultati e grafici i cui 
significati sono riassunti 
#nei remarks. 
 

Ho scritto le precedenti routines che sembrano funzionare, come si vede dai risultati,  considerando il periodogramma come una function, una specie di subroutine. Sarò costretto comunque a rimettere in discussione con calma altri programmi in R che contengono questa function tenendo conto dei cambiamenti!

 

CHI VOLESSE PUO’ VEDERE ANCHE GLI SCRIPTS DELLO STESSO AUTORE RELATIVI AL PERIODOGRAMMA E ALL’ANALISI DI FOURIER IN MATHEMATICA DI WOLFRAM VERS. 4.2, per fare un controllo dei risultati. Sono inseriti nelle appendici.

IL CONTROLLO  DEI PROGRAMMI IN R CHE SEGUONO E’ QUASI COMPLETATO

AD MAIORA

CHIUSA PARENTESI

________________________

period10001

6_CENNO A COMANDI IN R DI CALCOLO E ORGANIZZAZIONE DEI DATI

Filter, matrix e ts di R.
Discussione sui comandi di calcolo ed organizzazione sui dati. Commento sulle prime istruzioni di R (file di dati). Processi per automatizzare i “i conti”.

Si usa la funzione ts di R che riorganizza direttamente la serie originale (yt o as1)
in 12 colonne (mesi) e 5 righe (anni) per il calcolo poi con un for   delle medie di tutti i Gennaio, di Febbraio…

 

Discussione su filter

Applico direttamente la funzione Filter di R, sempre sulla serie originale (yt o as1), che, eliminando da essa (cioè da as1) la componente stagionale di ordine 12 + random, cambia contenuto in TREND + Ciclo + random? (divenendo la asf12).  Trovo poi la retta di regressione su asf12, i cui valori delle sue ordinate verranno tolti dalla serie originale; faccio il grafico di asf12 + retta di regr . Da controllare meglio. Smussando la yt, la asf12 è senza random? Vedere dopo gli script.

SEGUE IL COMMENTO SULLE  LE PRIME ISTRUZIONI DI R PER AUTOMATIZZARE I ‘CONTI’ DEL PROCESSO RIASSUNTO IN PRECEDENZA CHE ESPANDEREMO IN UN SECONDO TEMPO

I PRIMI INTERVENTI IN R

I primi passi nella schermata iniziale di R consistono nel caricare le Librerie suppletive di R necessarie a fornire i comandi, oltre a quelli di base, per gestire ed elaborare   i dati sperimentali. Con la funzione getwd() capisco dove ‘guarda’ R (cioè qual è la directory di lavoro) per cercare il file-dati da caricare e la funzione setwd (directory) permette di cambiare tale directory di lavoro. Fatta conoscere ad R la directory di lavoro, gli facciamo leggere il file-dati scelto per l’analisi (con il comando read.csv); nella fattispecie “As-Carlina1.csv”; la funzione file.show(“nome file.csv”) permette di visionare il contenuto del file che in generale è una matrice con righe e colonne è cioè un data.frame a cui si attribuisce un nome (per es., frame) e di cui è possibile conoscere le dimensioni col comando dim() o estrarre elementi. Le righe della matrice sono le osservazioni o casi; le colonne sono i campi o variabili. Con frame$variable si vuol dire di estrarre la variabile chiamata variable dal data.frame chiamato frame; frame[1,] significa prendere la prima riga, mentre frame[,3], prendere la terza colonna e così via. L’espressione summary(frame$variable) trova tutti i valori della variabile variabile contenuti nel data.frame chiamato frame. Così summary(frame[,3]), trova tutti i valori della colonna 3.

library (stats)

library(tseries)

library(lattice)

#library(graphics)

 

#getwd()

#setwd(“E:/R-2.12.2/bin/i386”)

# Se conosco dove è memorizzato il file con i dati da analizzare e la sua struttura

# utilizzo questi scripts iniziali

#as=read.csv(“As-Carlina.csv”)

#as1=as[,5]

#leggo la quinta colonna del data.frame: As-Carlina.csv dove c’è appunto yt

#as1=ts(as1) # considero as1 una serie storica

#ts.plot(as1) # plotto as1

Introdurremo invece direttamente la Serie yt o as1

as1= c(.033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043,

.066,.053,.082,.06,.08,.076,.056,.036,.05,.053,.056,.058,

.061,.063,.065,.068,.0815,.095,.079,.063,.069,.074,.08,.0765,.073,

.0695,.066,.093,.083,.073,.063,.074,.067,.06,.086,.08,.073,.067,

.089,.064,.087,.079,.07,.065,.06,.063)

7 – ECCO QUELLO CHE FAREMO CON R: ‘LETTURE’ SUI PROCESSI (‘CACCIA AI RESIDUI’ compresa)

ECCO QUELLO CHE FAREMO CON R

RIORGANIZZAZIONE DELLA SERIE STORICA MENSILE LUNGA CINQUE ANNI, As1, IN DODICI COLONNE (mesi)  E CINQUE RIGHE (anni) E BREVI LETTURE SUCCESSIVE

Il primo passo è riorganizzare la serie storica mensile della durata di 5 anni (5×12=60 mesi), in 12 colonne (mesi) e 5 righe (anni).

In ogni colonna ci sono 5 valori di ogni mese: nella prima, i 5 valori di gennaio, nella seconda, i 5 di febbraio e così via, Questo insieme costituisce il file as1.ts1. Per costruire as1.ts1 si può con R operare in almeno due modi. Una volta costituita la classificazione as1.ts1, si usa la funzione ts che permette poi tramite la subas, di meccanizzare con un for il calcolo delle dodici medie riferite ad ogni mese per i 5 anni (vedere dopo).

In sintesi con ts, che ha come argomenti: file, start e frequency, raggruppo i dati con i valori di ogni mese nella stessa colonna. Nella tabella appaiono il nome dei mesi su ogni colonna e il nome degli anni ad ogni riga; siamo così in grado di prendere i cinque dati di ogni mese (uno ogni dodici) per farne la media.

as1.ts1=ts(as1,start=1989,frequency=12)

Questa espressione fa anche la media di ogni colonna?

subas=as1.ts1[seq(1, length(as1), by=12)]

subas raccoglie i dati di gennaio per i 5 anni e ne fa la media(0.064); per ulteriori elaborazioni si può automatizzare con for.

Con for ottengo le 12 medie di ogni mese per 5 anni, mettendo un i al posto di 1 nell’argomento.

Guardiamo come.

mediamesi=c()

for(i in 1:12){mediamesi[i]=mean(as1.ts1[seq(i,length(as1),by=12)])}

ts.plot(mediamesi)

Se togliamo dal vettore mediamesi la media di as1, si ottiene una sorta di Effetto Stagionale mensile.

Mediamesi0=c()

Mediamesi0 =(mediamesi – mean(as1)) # da errore!

ts.plot(mediamesi0) # da errore! In effetti (vedere gli scripts al termine), non so perchè, sono necessarie variabili intermedie.

 

Vedremo dopo altri modi per il calcolo dell’Effetto Stagionale attraverso una Media Mobile e la funzione filter su as1, ambedue di ordine 12, modificando la stessa as1 o yt, in Mbt e asf12 di 12 termini più corte rispettivamente, contenenti ambedue almeno TREND lin.+ Ciclo (il random plausibilmente si cancellerebbe nel processo). La serie originale era pensata costituita da componente stagionale + TREND_ lin. + ciclo + random.

Calcolo la Media Mobile di ordine 12 su yt o as1; trovo la serie Mbt di 12 termini più corta, che è yt smussata della stagionalità, che serve a calcolare l’Effetto Stagionale, passando attraverso la sottrazione yt – Mbt , chiamata STRD (stagionalità più random:  Tabella 1, colonna 7, APPENDIX 3).

yt=as.vector(yt): n=length(yt); Mbt=c()

for(t in 7:n){Mbt[t] = (yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/12}

Mbt # di 12 termini più corta: 6 NA all’inizio e 6 NA alla fine, in tutto 48 dati (yt o as1 erano 60)

Mbt=Mbt[7:54]# elimino da Mbt gli NA; se i dati iniziali iniziavano da gennaio, Mbt inizia da un luglio e termina a un giugno

In alternativa applico il filter di ordine 12 su as1 o yt:

asf12=filter(yt, filter=rep(1/13,13)) # 12 o 13?

asf12

asf12=asf12[7:54] # elimino da asf12 gli NA

Le deboli differenze fra Mbt e asf12 è facile siano dovute alla Media Mobile manuale che è centrata.

Scorcio la as1 di 6 valori iniziali e finali per renderla lunga come Mbt e poi vi sottraggo Mbt:

STRD=as1[7:54] – Mbt # il primo valore di STRD corrisponde a luglio del primo anno.

Ciò significa: STRD= (ciclo+TREND+stagionalità+random) – (ciclo+TREND)=stagionalità+random; 60-12=48 termini.

Si calcola ora il Fattore Stagionale mensile (Tabella 1, colonna 8; 12 termini, APPENDIX 3) agendo con la funzione matrix su STRD e successivamente con colMeans: metto STRD (48 termini) sotto forma di matrice con dodici colonne (mesi) e 4 righe (anni)

stag = matrix(STRD, ncol=12, byrow=T)

Su questa matrice col comando colMeans posso trovare le 12 medie dei 4 valori, una per ogni mese, che metto in mediacol:

mediacol = colMeans(stag) # in mediacol rimangono i random?

Ordino le 12 medie ottenute, che iniziano da luglio del primo anno e terminano a giugno dell’anno successivo, da gennaio a dicembre:

mediacol=(mediacol[7:12],mediacol[1:6]) # Controllare se funziona!

mediacol # detto talora Fattore Stagionale

Copro poi i 5 anni ripetendo questi 12 valori:

ESAs = rep(mediacol,5) # Effetto stagionale di yt o as1

ESAs # serie lunga come yt o as1 originale

Dobbiamo ora togliere da yt o as1 l’Effetto Stagionale trovato per ottenere la serie iniziale destagionalizzata (stg, detta anche y1t o dst; Tabella 2, colonna 2) :

stg=c() #forse è meglio chiamala dst o y1t al posto di stg

dst=c() # dst o y1t in stg!

dst= yt–ESAs # TREND+ciclo_random; serie originale destagionalizzata (GRAF. N.4 a- CORR; b-PERIOD))

# Di fatto questa istruzione stranamente dava errore; forse è necessario introdurre variabili intermedie (vedere scripts relativi dopo). Controllare meglio!

# dst <- c(as1–ESAs) # TREND+ciclo_random #ancora da rifletterci!

dst  # è la serie originale destagionalizzata (in altre occasioni chiamata y1t). Di questa disegno il correlogramma: i dati sono autocorrelati; la statistica  DW , per K= 1,   N=60, rischio 0.05, cade a sinistra dell’intervallo dl-1.62 e si intravede la presenza di un TREND positivo (GRAF. N.4 a); dal periodogramma è sparito completamente il picco di frequenza 5 (periodo 60/5) dell’oscillazione stagionale (GRAF. N.4 b), presente invece  nel periodogramma della serie originale (GRAF. N.2 b) e nell’ESAs (GRAF. N.3 b).

y1t=dst

period0002

 

6-7 LA ‘CACCIA’ AI RESIDUI

Potremmo tentare di togliere da dst o y1t (TREND+ciclo_random) i random, provando a perequare con una Media Mobile 3*3 (pesata 1,2,3,2,1) per cui l’yt_smussato verrebbe a contenere ciclo+TREND che, tolto da dst o y1t, dovrei ottenere i random, se le ipotesi iniziali fossero giuste (vedere il testo di questi  scripts già in Blocco Note con  i risultati relativi, nel paragrafo prima delle Appendici (SECONDA PARTE). Alcuni ricercatori infatti propongono medie mobili a tre o 5 termini pesati 12321, per eliminare i random! PROVIAMO  invece il tentativo più classico che Segue: detrendizziamo linearmente la dst o y1t, sottoponendola ad una regressione lineare semplice (RLS)…

 

8 – INIZIO COPIA SCRIPTS DEL PROGRAMMA CENTRALE
Vari commenti possibili e riflessioni alternative

INIZIANO GLI SCRIPTS DEL PROGRAMMA RELATIVO A TUTTO IL PROCESSO DESCRITTO E DISCUSSO IN PRECEDENZA

Da copiare sul Blocco Note con copia/incolla e poi sulla consolle di R (o direttamente su R). In generale i programmi scritti in R o si fanno girare scrivendo una istruzione dietro l’altra , oppure, per es., si copiano gli  scripts sul Blocco Note od altro semplice programma di scrittura (anche quelli indirizzati ad R),  con copia/incolla e poi  sulla consolle di R.

Altro problema in R,  quando si copiano programmi pronti dal Blocco Note, è quello di gestire la visione dei diversi grafici, man mano che il programma gira. In questo caso è necessario che il programma controlli i grafici nel senso, per es., di far fermare il programma all’apparire del grafico nella finestra grafica, nella attesa della pressione di un tasto. Per questo esiste un semplice comando, da inserire, per es., all’inizio degli scripts, che ha la sintassi: par(ask=T).  Si può utilizzare in alternativa o insieme il comando par(mfrow=c(x,y) , che divide l’unica finestra grafica in x*y parti; x=2 e y=3, la finestra rimane divisa in 6 parti e può contenere 6 grafici e così via.

 

COMMENTO

Il seguente programma è stato utilizzato da prima nell’analisi della serie As originale, nel modo come era nato, cioè iniziando il lavoro con l’applicare la media mobile direttamente sulla serie originale, arrivando però ad una serie residuale  che può non rispettare i criteri richiesti (rivedremo i passaggi). Questo primo modo  è quello che per ora continua  a venire presentato e commentato.

Per osservare il percorso che parte invece, forse più giustamente, dalla serie detrendizzata (il trend in una serie  può  ‘disturbare’ il computo dell’Effetto Stagionale?), basta sostituire nel vettore as1, invece dei valori originali, i valori della serie detrendizzata, nel nostro caso per es. copiati dai programmi del Mathematica di Wolfram (Appendix 5) o dall’altro post  ‘Verso il periodogramma’, sempre dello stesso autore o… si rifaccia il conto. Basta togliere il cancelletto (#) all’as1 che riporta i valori della serie detrendizzata e ‘cancellettando’ invece i valori  dell’as1 che riporta  quelli della serie originale (e viceversa). I risultati ipoteticamente dovrebbero migliorare. Proviamo.

RESIDUI ANALISI SU As1 DETRENDIZZATO

Col tempo e la pazienza è possibile che riporti, in un link, il programma in pdf  che, in as1, ha i suoi valori detrendizzati, con più di una  decina di grafici relativi, con risultati e commenti! Vedere sopra la prima versione.

8-INIZIO COPIA PROGRAMMA

library(tseries)

library(lattice)

library(graphics)

as1= c(.033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043,.066,.053,

.082,.06,.08,.076,.056,.036,.05,.053,.056,.058,

.061,.063,.065,.068,.0815,.095,.079,.063,.069,.074,.08,.0765,.073,

.0695,.066,.093,.083,.073,.063,.074,.067,.06,.086,.08,.073,.067,

.089,.064,.087,.079,.07,.065,.06,.063)

# Per partire con la detrendizzazione, ad as1 sostituiamo i valori della stessa serie detrendizzata.

# Togliamo il cancelletto e mettiamo la nuova serie detrendizzata  qui e ‘cancellettiamo’ la precedente:

#as1 =c(-.018,.0089,-.0013,.0062,.0077,.0093,

#-.0012,-0.0187,-.0091,.00039,.0069,-.0085,

#-.0040,-.014,.0080,-.0054,.0231,.00064,

#.0202,.0157,-.0048,-.0252,-.0117,-.0092,

#-.0066,-.0051,-.0026,-.0011,.00048,.0030,

#.0160,.029,.013,-.0039,.0017,-.0092,

#.012,.0076,.0038,-.00018,-.0042,.0223,

#.012,.0014,-.0090,.0015,-.0060,-.0134,

#.0121,.0056,-.0018,-.0083,.0132,-.00122,

#.0102,.0018,-.0077,-.0131,-.0186,-.0161)

as1=ts(as1)

par(ask=T)

par(mfrow=c(1,2))

yt=c()

yt=as1

ts.plot(yt, main=”GRAF. N.2_yt_ SERIE CORRETTA”)

lines(yt,type=”l”)

acf(yt, main=”GRAF. N.2_a-yt_CORR_SERIE CORRETTA”)

#alfa=-pi/2 -> 270°; alfa=-1.175 rad (cioè -100°) -> 260°

 

#INIZIO FUNCTION

PRDGRAM<- function(y1,n1,m1) {

# VALORI DEL PARAMETRO ak

a0=c(); k=0; a0=0;

for(t in 1:n1){a0=a0+y1[t]*cos(2*pi*t*k/n1)}

a0

a0=a0*2/n1;a0=a0/2

a0

a=c();a[1:m1]=0;

for(k in 1:m1) {

for(t in 1:n1){

a[k]=a[k]+y1[t]*cos(2*pi*t*k/n1)}}

a=2*a/n1

# vALORI DEL PARAMETRO bk

b=c();b[1:m1]=0;b0=0;k=0

for(k in 1:m1) {

for(t in 1:n1){

b[k]=b[k]+y1[t]*sin(2*pi*t*k/n1)}}

a <- as.vector(a)

for(i in 1:m1){

if (abs(a[i]) < 1e-10) a[i]=0 else a[i]=a[i]}

a

for(i in 1:m1){

if (abs(b[i]) < 1e-10) b[i]=0 else b[i]=b[i]}

b=2*b/n1

b

# AMPIEZZE

#ro[1:m1]=0

ro <- sqrt(a^2 +b^2)

for(i in 1:m1){

if (abs(ro[i]) < 1e-10) ro[i]=0 else ro[i]=ro[i]}

# CALCOLO DELLA FASE DI OGNI ARMONICA

# RIPORTANDO IL VALORE AL QUADRANTE GIUSTO

f2=c()

f2[1:m1]=0

for(i in 1:m1){

f2[i] <- abs(a[i]/b[i])

f2[i] <- atan(f2[i])*180/pi}

f2 =as.vector(f2)

f2

#f2[1:m1]=0 un f2[1:m1] di troppo!

phi <- c()

for(i in 1:m1){

# f2 <- abs(a[i]/b[i]);

# f2 <- atan(f2)*180/pi;

if(b[i]>0 & a[i]>0) phi[i] = f2[i];

if(b[i]<0 & a[i]>0) phi[i] = 180-f2[i];

if(b[i]<0 & a[i]<0) phi[i] = 180+f2[i];

if(b[i]>0 & a[i]<0) phi[i] = 360-f2[i];

if(b[i]==0 & a[i]==0) phi[i] = 0;

if((b[i]<0 & b[i]>0) | a[i]==0) phi[i]=0;

if(b[i]==0 & a[i]>0) phi[i]=90;

if(b[i]==0 & a[i]<0) phi[i]=360-90

}

# PHI FASE ARMONICHE

phi=as.vector(phi)

phi

param_a <-a

param_b <-b

ampiezza <- ro

fase <- phi

# Qui, al termine della function si pone il valore di un’unica

# variabile che esce o, se escono più variabili, si usa

# un data.frame: data=data.frame(x1,x2,…).

# Ogni chiamata alla function permette di includere l’unica

# variabile o i data nel nome della chiamata:

# es. periodxx=nome.function(x1,x2,…)

data <-data.frame(a,b,ro, phi)

data

# questa matrice esce dalla function e viene ‘raccolta’ nella variabile nomexx (es.,periodxx)

}

#FINE FUNCTION

#Per richiamare la function:

#nomexx = PRDGRAM(Nome_var_vettore dati, numerosità del campione, numero di armoniche da cercare)

yt=as1

yx=as1

nx=length(yt)

#periodogramma yt

if (nx/2== nx%%2) mx=nx/2-1  else mx=(nx-1)/2 #da controllare se non sia necessario uno swap!

period_as1= PRDGRAM(yx, nx ,mx)

#par(mfrow=c(1,4)) 
#plot(a, xlab="Armoniche = N° osc. in n dati") 
#plot(b, xlab="Armoniche = N° osc. in n dati")

 period_as1 # tabella dei dati in uscita: ak, bk, ampiezze, fasi
# Con questa tabella si costruiscono le formule analitiche delle armoniche

period_as1$ro # vettore delle ampiezze

plot(period_as1$ro,type="l",main="GRAF. N.2; a-period_yt", 
xlab="Armoniche = N° oscill. in n dati", ylab="ampiezza")
 






As1_Corr_graf


period_su_As0001

 

par(mfrow=c(1,4))

plot(period_as1$a,ylab="Parametro a")
plot(period_as1$b,ylab="Parametro b") 
plot(period_as1$ro,type="l",main="PERIODOGRAMMA di as1", 
xlab="Armoniche = N° osc. in nx dati", ylab="ampiezza") 
plot(period_as1$phi,type="l", ylab="Fase")

#Per vedere i risultati trasferiti dalla consolle di R in pdf
#del precedente frammento di programma cliccare sotto:
As1_corr_R

par(mfrow=c(1,1)) 

as1.ts1=ts(as1,start=1989,frequency=12)
subas=as1.ts1[seq(1,length(as1),by=12)]

#-----------------------------------------------

# Gli scripts che riguardano il calcolo delle variabili vettoriali mediamesi e Mmesio per ora sono esclusi.

#mediamesi=c()

#for(i in 1:12){mediamesi[i]=mean(as1.ts1[seq(i,length(as1),by=12)])}

#ts.plot(mediamesi,main”mediamesi in 5 anni”)

#Mmesi0=c()

#a=mediamesi

#b=mean(as1)

#c=a-b

#Mmesi0=c () 12 valori medi meno la media serie originale; una specie di Effetto Stagionale

#Mmesi0=mediamesi – mean(as1)

#ts.plot(Mmesi0) # da controllare: Effetto Stagionale da confrontare con mediacol

#acf(Mmesi0, main=”CORR_Mmesi0″)

#Mmesi0 # da confrontare con mediacol

#—————————————————————————–

yt=as1

yt=as.vector(yt);  n=length(yt); Mbt=c()

for(t in 7:n){Mbt[t] = (yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+

yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/12}

#SI LAVORA ORA SU Mbt

Mbt #è quello che resta di as1, dopo la media mobile 12 (trend-ciclo_random)

Mbt=Mbt[7:54]# elimino da Mbt gli NA; Tabella N.1, colonna 6.

ts.plot(Mbt, main=”GRAF. N.4′; Mbt )

acf(Mbt, main=”GRAF. N.4′; acf_Mbt”)

#Periodogramma Mbt, serie più corta senza stagionalità

y3=c()

y3=Mbt

n3==length(y3)

if (n3/2== n3%%2) m3=n3/2-1  else m3=(n3-1)/2

#ifelse(nx%%2 > 0, m=(n-1)/2, m=n/2-1

period_Mbt=PRDGRAM(y3, n3 ,m3)

period_Mbt # tabella ak, bk,ro,phi

period_Mbt$ro #valori ampiezza di Mbt

ts.plot(period_Mbt$ro, main=”GRAF. N.4′; period_Mbt”)

 

# Filtro col comando filter la serie yt

asf12=filter(yt, filter=rep(1/13,13))

asf12

asf12=asf12[7:54] # elimino da asf12 gli NA

#Mbt  contiene l’as1 senza la stagionalità; in as1 però rimane quello

#che aveva ( trend-stagionalità-ciclo_random); se da as1, tolgo as1 senza la stagionalità,

#trovo la stagionalità e random (STRD) che trasformo in Effetto Stagionale eliminando

#una buona parte dei random.

FINE OPERAZIONI SU Mbt

#INIZIO CALCOLI CHE PORTANO ALL’EFFETTO STAGIONALE

STRD=as1[7:54]-Mbt # componente stagionale + random, serie più corta

STRD # da essa si estraggono gli Effetti Stagionali; TABELLA N.1, colonna 7:APPENDIX 3.

#Processo per costruire gli Effetti Stagionali attraverso STRD

stag = matrix(STRD, ncol=12, byrow=T) # variabile di passaggio a mediacol

mediacol = colMeans(stag) #in mediacol rimangono i random? o si perdono nella mediazione; 12 valori osc. annuale.

# in questo primo mediacol ottengo 12 valori a partire da luglio; TABELLA N.1, colonna 8; APPENDIX 3.

mediacol=c(mediacol[7:12], mediacol[1:6]) # qui ordino da gennaio a dicembre i 12

#valori dell’ EFFETTO STAGIONALE;

mediacol # è detto anche Fattore Stagionale; TABELLA N.1, colonna 8; APPENDIX 3.

#ts.plot(mediacol) # L’oscillazione annuale che copre 12 mesi (max in luglio)

ESAs = rep(mediacol,5) # l’Effetto Stagionale che ‘copre’  i 60 dati di yt o as1

ESAs #serie lunga come yt o as1 originale; TABELLA N.2, colonna 1; APPENDIX 3.

ts.plot(ESAs,main=”GRAF. N.3′; EFFETTO STAGIONALE”)

 

ESAs1 = rep(mediacol,2)

ts.plot(ESAs1,main=”GRAF. N.3; a-“EFFETTO STAGIONALE RLS”) #2 ascillazioni

acf(ESAs1, main=”GRAF. N.3′; b-CORR_EF. STAG. 2 ripet”)

#periodogramma ESAs1

yes=ESAs1

nes=length(ESAs1)

if (nes/2== nes%%2) mes=nes/2-1  else mes=(nes-1)/2

period_ESAs1=PRDGRAM(yes, nes, mes)

period_ESAs1

period_ESAs1$ro

plot(period_ESAs1$ro,type=”l”, main=”GRAF. N.3; b-Period_ro EFFETTO STAG.”)

 

dst=c() #attivo la serie destagionalizzata; dst o y1t ; TABELLA N.2, colonna 2; APPENDIX 3.

dst=as1-ESAs # da provare se funziona; destagionalizza

dst

#e=as1

#f=ESAs

#g=e-f

#dst=g

#Potrei smussare dst con una Media Mobile Pesata (3*3, cioè con  pesi 1,2,3,2,1) per tentare

#di eliminare la componente casuale

#Si otterrebbe una serie (y1t) contenente CICLO+TREND, che se la tolgo dalla serie destagionalizzata

#dst precedente dovrei ottenere il #RESIDUO.

yd=dst

nd=length(dst)

if (nd/2== nd%%2) md=nd/2-1  else md=(nd-1)/2

period_dst=PRDGRAM(yd, nd, md)

period_dst

period_dst$ro

plot(period_dst$ro,type=”l”,main=”GRAF. N.4: b-Period. dst o y1t”)

#PROVIAMO INVECE A TOGLIERE IL TREND DALLA dst o y1t

plot(dst,type=”l”, main=”yt-destagionalizzata”) # la y1t o dst= yt destagionalizzata= ciclo+TREND +random (GRAF. N.4′)

acf(dst, main=”GRAF. N4; a-CORR-y1t o dst”)

# Se elimino il TREND da dst ottengo CLRD e posso controllare con CORR se

# ciò che resta è da considerare residuo. yt-ESAs-TREND = CLRD

# CLRD =yt-TREND- ESAs

#Calcolo il trend di dst per toglierlo da yt-ESAs o da y1t ed ottenere CLRD

t=seq(1:60)

fitdst=lm(dst~t)

abline(lm(dst~t))

summary(fitdst)

resid(fitdst)

p=predict(fitdst,data.frame(t=c(1,60)))

CLRD=c()

CLRD=dst-p

CLRD

CLRD=yt-ESAs-p

n1=length(p)

ts.plot(CLRD, main=”GRAF. N.5-RESIDUI” )

acf(CLRD, main=”GRAF. N.5; a-CORR_CLRD”)

#periodogramma di CLRD

yr=CLRD

nr=length(yr)

if (nr/2== nr%%2) mr=nr/2-1  else mr=(nr-1)/2

period_clrd=PRDGRAM(yr, nr ,mr)

period_clrd

period_clrd$ro

plot(period_clrd$ro,type=”l”,main=”GRAF. N.5; b-Period. CLRD”)

#da controllare ancora!

#FINE COPIA PROGRAMMA  da trasferire in Blocco Note o direttamente sulla consolle di R

PER VEDERE SCRIPS E COMMENTI PRECEDENTI + RESULT  IN pdf CLICCARE SOTTO:

ANALI SU As1 DETRENDIZZATO

BLOCCO_NOTE_PERCORSO_PERIOD0

BLOCCO_NOTE_PERCORSO_PERIOD

COMMENTO

Sembra che in questo processo CLRD (residui) non siano random e siano correlati (da provare altri tests. Proviamo però a fare altre misure di controllo. Se è così percorriamo altre vie già accennate. Possiamo partire col detrendizzare la serie originale as1, rendendola nelle previsioni stazionaria, e procedere con gli stessi scripts già usati.

Se ai dati originali di as1  sostituiamo i  dati originali senza però il trend rettilineo (serie originale detrendizzata, nelle previsioni resa stazionaria), possiamo vedere che cosa accade. In effetti sembrerebbe che, se invece partiamo coll’applicare  una media mobile di ordine 12 su una serie non stazionaria, si possa arrivare a questo risultato.

Se si parte con una detrendizzazione (serie stazionaria) e poi si applica la media mobile per trovare gli Effetti Stagionali, che togliamo dalla serie originale, e si procede con successiva detrendizzazione su serie_originale- Eff. Stag., si prevede un aumento dell’ R-quadro e forse un risultato più idoneo.

Si fa prima una regressione sulla serie di partenza; attraverso una media mobile si cercano gli Effetti Stagionali che togliamo dalla serie originale (la non stazionarità può disturbare gli effetti stagionali), ottenendo la serie originale destagionalizzata;  si fa infine una seconda regressione su questa differenza, cioè sulla serie destagionalizzata, che può  contenere appunto TREND + CICLO_RANDOM, ricavando poi il CICLO_RANDOM (da verificare).

Altro percorso: analisi dei dati trimestrali della stessa serie as1.

9 – PRIMA PARTE IN SINTESI

PRIMA PARTE IN SINTESI

LA SERIE PEREQUATA Mbt, L’EFFETTO STAGIONALE ESAs, LA SERIE DESTAGIONALIZZATA y1t (dst), LA y1t SMUSSATA: ciclo+TREND (y1ts),

LA COMPONENTE CASUALE O RESIDUI

 IL CORRELOGRAMMA, IL TEST DI DURBIN WATSON  e di LINMUDHOLKAR

Dopo aver eliminato la componente stagionale (ESAs : APPENDIX3, TABELLA N.2, col.1) dalla serie originale yt  (APPENDIX3,  TABELLA N.1, col.5) sottraendo yt – ESAs, si ottiene la serie destagionalizzata (dst ovvero y1t:  APPENDIX3, TABELLA N.2, col.2). In questa serie sanno rimasti gli eventuali ciclo, TREND e la componente random. Sottopongo quest’ultima al programma CORR : i dati sono autocorrelati positivamente (la statistica di Durbin Watson , per k= 1, N=60 e rischio 0.05, cade a sinistra dell’intervallo dl-du (1.55-1.62) e si nota la presenza un TREND positivo (GRAF. N.4 a); dal periodogramma è completamente scomparso il picco di frequenza 5 (periodo 60/5) dell’oscillazione stagionale (GRAF. N.4 b), presente invece nel periodogramma della serie originale (GRAF. N.2 b) e nell’ESAs (GRAF. N.3 b). Leggere Appendice 1.

Smussiamo la serie y1t o dst con una media mobile pesata 3*3  (1,2,3,2,1), per eliminare la componente casuale. Si ottiene così la serie y1ts (CLTR : APPENDIX3, TABELLA N.2, col.3) che potrebbe contenere nelle previsioni ciclo e TREND (CLTR). Sottraendo da y1t o dst (ciclo+TREND+Random) la serie y1ts che potrebbe contenere ciclo+TREND si dovrebbe ottenere la componente casuale o serie random. Testando tale serie col programma CORR, risulta che essa è rumore di fondo (white noise), avvalorando il processo usato fino a questa fase. Infatti la DW, per k=1, n=60 e alfa =0.05, ha valore 2.57 (vedere tabella Appendice 2) per cui esce dall’intervallo ricavato dalle tabelle dl-du (1.55-1-62): assenza di correlazione interna. la statistica di LIN-MUDHOLKAR, per la gaussiana, per alfa=0.05 e r=+/- 0.403 ricavato dalle tabelle, ha il valore -0.0416, cioè cade all’interno dell’intervallo di r, per cui non posso rifiutare l’ipotesi nulla: la distribuzione dei residui così calcolati è da considerarsi gaussiana. Forse è proprio l’effetto di non aver esplicitata la serie CLTR  con il calcolo del TREND a favorire la compatibilità dei residui alle ipotesi iniziali.

 

MODELLO DI REGRESSIONE LINEARE SEMPLICE  (RLS) E TESTS RELATIVI.

ADEGUAMENTO DEL MODELLO DI REGRESSIONE ALLA POPOLAZIONE. COEFFICIENTI DELLA RETTA ED R-q

TEST SU R-q E LA F DI FISHER , TESTS SUI COEFFICIENTI DELLA RETTA, INTERVALLO DI CONFIDENZA.

RESIDUO DELLA REGRESSIONE E L’AFFIDABILITA’.

Applichiamo invece  a y1t o y1ts (APPENDIX3, TABELLA N.2, col.2;  TABELLA N.2, col.3) )  un modello di regressione per separare il TREND dai loro contenuti.  Proviamo una regressione lineare con la sola variabile, il tempo, misurato in mesi (un solo regressore, k1=1 nelle tabelle DW), senza preoccuparci per ora se tale modello sia idoneo. Lo controlleremo dall’analisi dei residui. Se sono rispettate le assunzioni di linearità, una buona misura dell’adeguamento del modello lineare ai dati è il Coefficiente di Determinazione R-quadro. La sua radice quadrata R è il Coefficiente di correlazione di Pearson detto anche Multiple-R. Se R-q è 1, significa che tutte le osservazioni cadono sulla retta di regressione; se  zero, nessuna associazione lineare fra le variabili, anche se può esserci una relazione non lineare. R-q può così essere interpretato come la proporzione della variazione di y ‘spiegata’ dal modello , come precisato in altre occasioni. Su y1t o su y1ts, si opera con una regressione lineare calcolando bo e b1 ed ottenendo in ambedue i casi, come era prevedibile, differendo le due serie per la sola componente casuale, la stessa retta di regressione seguente (APPENDIX3, TABELLA N.3, col.8 e  APPENDIX3, TABELLA N.4, col.3 per i valori previsti):

y_predetto = TREND = TREND’ = 0.051 + 0.00005*t

Vedere  APPENDIX3, TABELLA N.3, col.5, per i risultati intermedi al fine del cacolo dei coefficienti della retta.

Nel nostro caso  R-q = 0.44, cioè il modella spiega il 44% della variazione complessiva della variabile dipendente. Per controllare l’ipotesi  nulla che nella popolazione non esista relazione lineare (R-q_pop.=0), si procede con l’analisi della varianza. Per tutti i particolari dei ‘conti’ che seguono vedere, per es.,  il Post ‘Un percorso verso il periodogramma’ su questo stesso BLOG. Seguendo le indicazioni riportate nel paragrafo relativo a questo argomento nel Post  su nominato, si ottiene la seguente tabella:

                                                       GL          SOMMA DEI QUADRATI          MEAN SQUARE

Variazione di regressione       1                              0.o0435                                       0.00435

Variazione residuale                58                           0.00559                                       0.000096

TOT                                                                               0.00994

da cui: Somma quadrati reg./Somma quadrati tot = 0.44, cioè R-quadro.

La statistica  F di Fisher che permette di saggiare l’ipotesi nulla: R-quadro pop.=0, è 0.00435/0.000096 = 45.31, da cui, riportata sulle tavole con 1 e 58 gradi di libertà (GL), si ricava una significanza per F minore di 0.00001, per cui si respinge l’ipotesi nulla e nella popolazione esisterà con alta probabilità una relazione lineare.

Procedendo ancora a prove incrociate si può testare l’ipotesi che  b1_pop. =0; si calcola la statistica T per b1: pendenza/errore standard_pend, ottenendo ERb1=7.31*10^-5 e poichè b1=0.000492, risulta T=6.73, che dalle tabelle relative per 58 gradi di libertà (GL=N-2) si ha una significanza per T di 0.0000..<<0.05, per cui si respinge l’ipotesi nulla che la pendenza della popolazione sia zero (quindi esiste dipendenza lineare).

Procedendo, nell’intervallo di confidenza al 95% per la pendenza non potrà allora il valore zero. Infatti calcolando ESb1 come suggerito da altri interventi (0.000073), l’intervallo di confidenza al 95% per beta1 risulta (con 58 GL):

b1-1.96*ESb1 <=  beta1 <=  b1+ 1.96*ESb1

0.00492-0.00014 <= beta1 <= 0.00492+0.00014

0.00035 <= beta1 <=0.00063

Si vede chiaramente come i vari tests, se affidabili, confermano la presenza di un trend lineare nei dati.

Togliendo da y1t la serie del trend, si otterrà la serie CLRD ( APPENDIX3, TABELLA N.4, col.4) con l’eventuale ciclo + la componente casuale (random) I residui della regressione, per il modo con cui abbiamo proceduto, sono proprio i valori della serie CLRD. E’ prevedibile che questa serie, se davvero includerà una componente ciclica significativa,non risulterà rispetterà almeno qualche condizione fra quelle ipotizzate sui residui (indipendenza, varianza costante…). procederemo ad investigare questa serie sui residui. Applicando ad essi il programma CORR, otteniamo il grafico, GRAF. N.5 a) correlogramma) e 5b (periodogramma), il test per l’indipendenza di Durbin Watson e quello per la normalità di Lin Mudholkar. Il valore di DW è risultato 1.378, che (N=60, K’=1 e alfa =0.05) esce a sinistra dell’intervallo 1.55-1.62 e quindi l’autocorrelazione è positiva, mentre il test per la gaussiana (rischio 0.05, N=60, r=+/-0.403, fornisce rc=-0.0298, cioè all’inteno dell’intervallo, per cui non posso rifiutare l’ipotesi nulla (la serie ha distribuzione gaussiana). Graficando i residui standardizzati con la variabile pred pure standardizzata, si ottiene il   GRAF. N. 6 a dove non appaiono patterns evidenti. Dal GRAF. N.6 b invece, ottenuto riportando i residui per ogni unità di tempo, si evidenzia una qualche variazione della varianza dei residui (eteroscedasticità, variazione a clessidra). Allora i tests che fanno riferimento al comportamento della popolazione universo (in particolare gli F-tests) possono non essere affidabili e quindi incerto il modello di regressione usata.

Al termine dell’analisi con un modello di regressione lineare semplice, tenteremo ulteriori approfondimenti alla ricerca di un maggiore R-quadro, ma specialmente di una maggior concordanza dei residui alle condizioni iniziali (linearità, normalità, indipendenza, omoscedasticità).

 

 

 

 

stat_period_corr0001


 stat_period_corr0002

i

stat_reg_mlr_blog0001

10 – SECONDA PARTE IN SINTESI

SECONDA PARTE IN SINTESI: UN ALTRO TENTATIVO SULLA CACCIA AI RESIDUI (senza passare attraverso una regresssione)

SCRIPTS IN BLOCCO NOTE:  DA COPIARE DIRETTAMENTE SULLA CONSOLLE DI R

 

# Intanto trascriviamo nel vettore yt i 60 dati della conc. As da cui partire. Impariamo poi a calcolare con R gli altri 5 vettori dati che faranno parte dell'analisi della nostra serie
# reale e quindi della nostra esercitazione. Calcoliamo come primo vettore Mt (media mobile di ordine  12 su yt.

yt=c(.033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043,.066,.053,.082,.06,.08,.076,.056,.036,.05,
.053,
.056,.058,.061,.063,.065,.068,.0815,.095,.079,.063,.069,.074,.08,
.0765,.073,.0695,.066,.093,.083,.073,.063,.074,.067,.06,.086,.08,.073,.067,.089,.064,.087,.079,.07,.065,.06,.063)

t=1

#Come primo passo grafichiamo i dati e osserviamo se ci sono regolarità all'interno (trend, oscillazioni), precisiamo le ipotesi con un correlogramma ed un periodogramma, I dati sono mensili: Ipotizziamo comunque una oscillazione di periodo 12.

# Calcoliamo, come primo vettore, Mt (media mobile centrata e pesata di ordine 12 su yt).

yt=as.vector(yt) ; n=length(yt); Mt=c()
for(t in 7:n){Mt[t] = (yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+
yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/12}
Mtc=Mt[7:54]

mt=filter(yt,filter=rep(1/13,13))
# calcolo della Mm col comando filter di R: confrontare i due risultati
mt #OK

# in Mt ci sono i 48 (60-12) dati Media mobile di yt, da cui costruisco i 12 Fattori Stagionali (FStag) 
facendo la media dei 4 gennaio, dei 4 febbraio ecc. a partire da luglio, perchè Mt iniziava con luglio.
FSTag0 = matrix(Mtc, ncol=12, byrow=T)
# matrice di 4 righe (valori dei 12 mesi dei 4 anni) e 12 colonne con in ognuna le 4 conc. dei mesi dello stesso nome a partire da un luglio.
FStag1=colMeans(FSTag0)
#  in FStag1 trovo le 12 medie dei 4 mesi dello stesso nome (inizio luglio, fine giugno)
FStag=c(FStag1[7:12], FStag1[1:6]) # da controllare! Ordino da gennaio. OK
ESAs=rep(FStag,5) # EFFETTO STAGIONALE As
ESAs # 60 dati
Yt1=yt-ESAs # Ciclo+Trend+Random
Yt1 # 60 dati
Yt1c=Yt1[3:58]
Yt1s=c()
for(i in 1:60){Yt1s[i]=(Yt1[i-2]+2*Yt1[i-1]+3*Yt1[i]+2*Yt1[i+1]+
Yt1[i+2])/9}
Yt1s=as.vector(Yt1s) # smusso Yt1 con Mm 3*3

ns=length(Yt1s) # più corto di 4 elementi
Yt1s # yt1 senza random; cioè Ciclo+Trend

par(ask=T)

Yt1s=Yt1s[3:(ns-2)]

RD=Yt1c-Yt1s # forse si tratta solo di random: il Ciclo?

#Riportiamo in una tabella 1 5 vettori dell'analisi su yt

#data <- data.frame(t,yt,ESAs,Yt1,RD)

# Facciamo i 5 correlogrammi dei vettori trovati: yt, ESAs, Yt1, Yt1s, RD
coyt=acf(yt)
coyt
coESAs=acf(ESAs)
coESAs
coYt1=acf(Yt1)
coYt1s=acf(Yt1s)
coYt1s
coRD=acf(RD)
coRD
# Interessante abbinare il correlogramma con il periodogramma e da controllare i correlogrammi con il programmino scritto dall'autore


RISULTATI DEL PROGRAMMA PRECEDENTE (come si vede gira senza errori!)


> # Interessante abbinare il correlogramma con il periodogramma.
> # Intanto trascriviamo nel vettore yt i 60 dati della conc. As da cui partire. Impariamo poi a calcolare con R gli altri 5 vettori dati che faranno parte dell'analisi della nostra serie
 # reale e quindi della nostra esercitazione. Calcoliamo come primo vettore Mt (media mobile di ordine  12 su yt.

> 
> yt=c(.033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043,.066,.053,.082,.06,.08,.076,.056,.036,.05,
+ .053,
+ .056,.058,.061,.063,.065,.068,.0815,.095,.079,.063,.069,.074,.08,
+ .0765,.073,.0695,.066,.093,.083,.073,.063,.074,.067,.06,.086,.08,.073,.067,.089,.064,.087,.079,.07,.065,.06,.063)
> 
> t=1
> 
> #Come primo passo grafichiamo i dati e osserviamo se ci sono regolarità all'interno (trend, oscillazioni), precisiamo le ipotesi con un correlogramma ed un periodogramma, I dati sono mensili: Ipotizziamo comunque una oscillazione di periodo 12.
> 
> # Calcoliamo, come primo vettore, Mt (media mobile centrata e pesata di ordine 12 su yt).
> 
> yt=as.vector(yt) ; n=length(yt); Mt=c()
> for(t in 7:n){Mt[t] = (yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+
+ yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/12}
> Mtc=Mt[7:54]
> 
> mt=filter(yt,filter=rep(1/13,13)) # 13 o 12?
> # calcolo della Mm col comando filter di R: confrontare i due risultati
> mt #OK
Time Series:
Start = 1 
End = 60 
Frequency = 1 
 [1]         NA         NA         NA         NA         NA         NA
 [7] 0.05115385 0.05192308 0.05369231 0.05384615 0.05561538 0.05553846
[13] 0.05684615 0.05861538 0.06015385 0.05938462 0.05892308 0.05815385
[19] 0.05876923 0.05915385 0.06053846 0.06030769 0.06123077 0.06015385
[25] 0.06180769 0.06296154 0.06319231 0.06373077 0.06626923 0.06811538
[31] 0.07019231 0.07176923 0.07292308 0.07357692 0.07380769 0.07596154
[37] 0.07711538 0.07646154 0.07400000 0.07361538 0.07392308 0.07323077
[43] 0.07415385 0.07415385 0.07388462 0.07342308 0.07492308 0.07476923
[49] 0.07430769 0.07400000 0.07376923 0.07392308 0.07284615 0.07253846
[55]         NA         NA         NA         NA         NA         NA
> 
> # in Mt ci sono i 48 (60-12) dati Media mobile di yt, da cui costruisco i 12 Fattori Stagionali (FStag) facendo la media dei 4 gennaio, dei 4 febbraio ecc. a partire da luglio, perchè Mt iniziava con luglio.
> FSTag0=matrix(Mtc, ncol=12, byrow=T)
> # matrice di 4 righe (valori dei 12 mesi dei 4 anni) e 12 colonne con in ognuna le 4 conc. dei mesi dello stesso nome a partire da un luglio.
> FStag1=colMeans(FSTag0)
> #  in FStag1 trovo le 12 medie dei 4 mesi dello stesso nome (inizio luglio, fine giugno)
> FStag=c(FStag[7:12], FStag1[1:6]) # da controllare! Ordino da gennaio. OK
> ESAs=rep(FStag,5) # EFFETTO STAGIONALE As
> ESAs # 60 dati
 [1] 0.06147115 0.06221154 0.06285577 0.06317308 0.06323077 0.06334615
 [7] 0.05878846 0.05965385 0.06022115 0.06050962 0.06085577 0.06113462
[13] 0.06147115 0.06221154 0.06285577 0.06317308 0.06323077 0.06334615
[19] 0.05878846 0.05965385 0.06022115 0.06050962 0.06085577 0.06113462
[25] 0.06147115 0.06221154 0.06285577 0.06317308 0.06323077 0.06334615
[31] 0.05878846 0.05965385 0.06022115 0.06050962 0.06085577 0.06113462
[37] 0.06147115 0.06221154 0.06285577 0.06317308 0.06323077 0.06334615
[43] 0.05878846 0.05965385 0.06022115 0.06050962 0.06085577 0.06113462
[49] 0.06147115 0.06221154 0.06285577 0.06317308 0.06323077 0.06334615
[55] 0.05878846 0.05965385 0.06022115 0.06050962 0.06085577 0.06113462
> Yt1=yt-ESAs # Ciclo+Trend+Random
> Yt1 # 60 dati
 [1] -0.0284711538 -0.0192115385 -0.0118557692 -0.0041730769 -0.0022307692
 [6] -0.0003461538 -0.0057884615 -0.0236538462 -0.0142211538 -0.0045096154
[11]  0.0021442308 -0.0131346154 -0.0084711538 -0.0192115385  0.0031442308
[16] -0.0101730769  0.0187692308 -0.0033461538  0.0212115385  0.0163461538
[21] -0.0042211538 -0.0245096154 -0.0108557692 -0.0081346154 -0.0054711538
[26] -0.0042115385 -0.0018557692 -0.0001730769  0.0017692308  0.0046538462
[31]  0.0227115385  0.0353461538  0.0187788462  0.0024903846  0.0081442308
[36]  0.0128653846  0.0185288462  0.0142884615  0.0101442308  0.0063269231
[41]  0.0027692308  0.0296538462  0.0242115385  0.0133461538  0.0027788462
[46]  0.0134903846  0.0061442308 -0.0011346154  0.0245288462  0.0177884615
[51]  0.0101442308  0.0038269231  0.0257692308  0.0006538462  0.0282115385
[56]  0.0193461538  0.0097788462  0.0044903846 -0.0008557692  0.0018653846
> Yt1c=Yt1[3:58]
> Yt1s=c()
> for(i in 1:60){Yt1s[i]=(Yt1[i-2]+2*Yt1[i-1]+3*Yt1[i]+2*Yt1[i+1]+
+ Yt1[i+2])/9}
> Yt1s=as.vector(Yt1s) # smusso Yt1 con Mm 3*3
> 
> ns=length(Yt1s) # più corto di 4 elementi
> Yt1s # yt1 senza random; cioè Ciclo+Trend
 [1]            NA            NA -1.255983e-02 -6.694444e-03 -3.708333e-03
 [6] -4.989316e-03 -9.090812e-03 -1.287073e-02 -1.140385e-02 -8.274573e-03
[11] -5.727564e-03 -8.419872e-03 -9.424145e-03 -1.017735e-02 -4.337607e-03
[16] -1.027778e-03  5.958333e-03  8.455128e-03  1.157585e-02  6.129274e-03
[21] -2.070513e-03 -1.060791e-02 -1.194979e-02 -9.530983e-03 -5.979701e-03
[26] -3.955128e-03 -2.004274e-03 -2.777778e-05  3.902778e-03  1.089957e-02
[31]  1.874252e-02  2.179594e-02  1.809615e-02  1.216987e-02  1.027244e-02
[36]  1.208013e-02  1.424252e-02  1.326709e-02  1.032906e-02  9.861111e-03
[41]  1.273611e-02  1.806624e-02  1.824252e-02  1.524038e-02  1.026282e-02
[46]  7.836538e-03  7.827991e-03  9.913462e-03  1.368697e-02  1.393376e-02
[51]  1.377350e-02  1.130556e-02  1.384722e-02  1.478846e-02  1.779808e-02
[56]  1.546261e-02  1.159615e-02  5.836538e-03            NA            NA
> 
> par(ask=T)
> 
> Yt1s=Yt1s[3:(ns-2)]
> 
> RD=Yt1c-Yt1s # forse si tratta solo di random: il Ciclo?
> 
> #Riportiamo in una tabella 1 5 vettori dell'analisi su yt
> 
> #data <- data.frame(t,yt,ESAs,Yt1,RD)
> 
> # Facciamo i 5 correlogrammi dei vettori trovati: yt, ESAs, Yt1, Yt1s, RD
> coyt=acf(yt)
Aspetto per confermare cambio pagina...
> coyt


Autocorrelations of series ‘yt’, by lag

     0      1      2      3      4      5      6      7      8      9     10 
 1.000  0.541  0.395  0.223  0.302  0.221  0.330  0.281  0.150  0.102  0.150 
    11     12     13     14     15     16     17 
 0.248  0.255  0.308  0.197  0.099 -0.006  0.042 
> coESAs=acf(ESAs)
Aspetto per confermare cambio pagina...
> coESAs

Autocorrelations of series ‘ESAs’, by lag

     0      1      2      3      4      5      6      7      8      9     10 
 1.000  0.542  0.187 -0.111 -0.349 -0.460 -0.455 -0.423 -0.309 -0.102  0.146 
    11     12     13     14     15     16     17 
 0.433  0.800  0.434  0.150 -0.088 -0.276 -0.362 
> coYt1=acf(Yt1)
Aspetto per confermare cambio pagina...
> coYt1s=acf(Yt1s)
Aspetto per confermare cambio pagina...
> coYt1s

Autocorrelations of series ‘Yt1s’, by lag

    0     1     2     3     4     5     6     7     8     9    10    11    12 
1.000 0.908 0.757 0.610 0.519 0.457 0.396 0.326 0.271 0.256 0.277 0.317 0.335 
   13    14    15    16    17 
0.321 0.263 0.198 0.145 0.123 
> coRD=acf(RD)
Aspetto per confermare cambio pagina...
> coRD

Autocorrelations of series ‘RD’, by lag

     0      1      2      3      4      5      6      7      8      9     10 
 1.000 -0.308 -0.166 -0.187  0.222 -0.198  0.195  0.066 -0.089 -0.097  0.014 
    11     12     13     14     15     16     17 
 0.004 -0.029  0.147  0.043 -0.114 -0.071  0.046 
> # Interessante abbinare il correlogramma con il periodogramma: da fare.

L’EPILOGO

SEGUONO ULTERIORI APPROFONDIMENTI: 

APPLICAZIONE DI UNA REGRESSIONE LINEARE MULTIPLA (RLM) OPPORTUNA (variabili "dummy").
 
COME CALCOLARE LA F DI FISHER NELLE REGRESSIONI RLM.
 
COME CALCOLARE L'ERRORE STANDARD SUI COEFFICIENTI DI REGRESSIONE NELLA RML
 
COME SI APPLICA UNA REGRESSIONE LINEARE MULTIPLA PESATA (RLMP)
 
CHI VOLESSE ESERCITARSI SU ESEMPI RELATIVI AL CALCOLO MATRICIALE APPLICATO ALL'ANALISI DI DATI    SPERIMENTALI CERCARE IN QUESTO SITO "TIPS DI SCIENZA" (in particolare sui "conti" relativi alla   regressione lineare multipla (MLR). 
LO SCRITTO CHE SEGUE E' L'ULTIMA TRANCE DELL'ARTICOLO ORIGINALE CHE RIGUARDA GLI ULTERIORI        APPROFONDIMENTI,ELENCATI SOPRA, SCRITTO ANCORA DALLO SCRIVENTE, RIVISITATO E INTERPRETATO CON R   IN QUESTO POST. I RIFERIMENTI COME 1.1.2.2 ECC. RIGUARDANO RIMANDI A SUOI PARAGRAFI SPECIFICI. DATA LA NATURA A 'ZIBALDONE LEOPARDIANO DISPERSO' DI QUESTO LAVORO A GETTO ROBINSONIANO CI PROPONIAMO DI INSERIRE LA SECONDA PARTE DELL'ORIGINALE PRIMA DELLE APPENDICI. DOVREMMO SCANNERIZZARLO MEGLIO! 

11-L'EPILOGO

stat13_ridot

reg_MLR_blog0003

14reg_MLR_blog0004

15

stat10001

stati0002

reg_MLR_blog0007

reg_MLR_blog0008

stat_reg_mlr_blog0001

i

durbin watson0008

i

GRAF. N.9

durbin watson0009

i

durbin watson0010

12 -APPENDICE1

APPENDICE1

 

Il correlogramma ed il test di Durbin-Watson. ([3], 949-953)

Ammettiamo che il lettore conosca il Coefficiente di Correlazione lineare di Pearson, ovvero date N paia di osservazioni su due variabili X  e Y, tale Coefficiente di Correlazione  fra esse e dato:

r =  Σi(Xi- Xm)*(Yi – Ym)/SQR[ Σi(Xi-Xm)^2 * Σi(Yi – Ym)]

Quest’idea viene trasferita alle serie storiche per vedere se osservazioni successive sono correlate.

Date N osservazioni X1, X2,………Xn , in una serie storica discreta possiamo considerare N-1 paia di osservazioni (X1,X2), (X2,X3), . . . ,(X(n-1),Xn), le cui prime osservazioni di ogni paio costituiscono la prima variabile e le seconde, la seconda variabile. Se si applica la formula precedente, dove Xi sarebbe Xt e Yi sarebbe Y(t+1), mentre Xm sarebbe la media della prima variabile (da t=1 a t=N-1) e Ym sarebbe la media della seconda variabile (da t=2 a t=N,  in ambedue i casi il numero degli elementi sarebbe N-1. Si otterrebbe una formula complessa con due medie diverse che vengono invece calcolate ambedue sulla serie originaria di numerosita N. Si usa cosi la formula approssimata scritta sotto, estesa al caso in cui si voglia trovare la correlazione tra serie di osservazioni a distanza H fra loro (slittate di h termini o di lag h)

I coefficienti di auto-correlazione rh , dove h=0,1,2…q e q è minore ad uguale a (N-2)/2, sono coefficienti di correlazione, calcolati per ogni valore di h, che misurano la concordanza o la discordanza tra i valori di una serie storica e quelli della stessa però slittati di h unità di tempo (lag h), consentendo di analizzare la sua struttura interna, ossia i legami fra i termini della stessa ([8] 18-20).

 rh = Σi[(y(t)-ym)(y(t+h)-ym)]/[(n-h)*Σj(y(t)-ym)^2/n)] dove i va da t=1…n-h e j va da t=1 … n

in alcuni testi viene abolito il fattore n/(n-h).

Tale formula presenta la semplificazione di poter   utilizzare una media unica per le Yt (quella dei dati originali), presupponendo una situazione stazionaria ([8] pag.19 e [2] pag.133). In particolare r0 = 1 (lag h =0, nessun slittamento) e gli altri rh assumono valori fra +1 (completa concordanza) e -1 (totale discordanza). Il correlogramma è la rappresentazione grafica dei coefficienti di auto-correlazione in funzione degli slittamenti (lags h) e permette di vedere se la serie storica possiede qualche regolarità interna.

CENNI DI LETTURA DEI CORRELOGRAMMI

-I coeff. di autocorr. di dati random hanno distribuzione campionaria che può essere approssimata da una curva gaussiana con media zero ed errore standard 1//N. Questo significa che il 95% di tutti i coeff. di autocorr. , calcolati da tutti i possibili campioni estratti, dovrebbero giacere entro un range specificato da: zero +/- 1.96 errori standard. I dati cioè della serie saranno da considerarsi random se questi coefficienti saranno entro i limiti:

 

-1.96 (1/√n)≤ rh  ≤ +1.96 (1/√n);       la fascia dell’errore:   +/- 2/√n

 

Per l’interpretazione dei correlogrammi vedere ([8] 20-25) da cui ricaviamo le seguenti informazioni.

 

 

– Una serie storica completamente casuale, cioè i cui successivi valori sono da considerarsi tutti indipendenti fra loro (non correlati), tutti i valori di rh  (eccetto r0 che è sempre +1, correlazione della serie con se stessa) oscilleranno accidentalmente intorno allo zero entro la fascia dell’errore. Se l’idea iniziale era questa in effetti  5 su 100 valori di rh potrebbero superare la fascia dell’errore e se plotto il correlogramma, 19 su 20 valori di rh potrebbero cadere all’interno della fascia, ma ci si potrebbe aspettare che uno possa eesere significativo sulla media. Insomma anche se la serie è casuale, ogni tanto verso lag più elevati potrebbero apparire picchi significativi. Se abbiamo a che fare con un numero elevato di coefficienti, potrebbero apparire risultati non aspettati. Questo rende il correlogramma uno strumento di investigazione incerto.

 

– I coeff. di autocorr. per i dati stazionari (assenza di TREND) vanno velocemente a zero dopo il 4° o 5° lag di tempo e  sono significativamente diversi da zero per i primi lag. Anche su correlogrammi,  ai lags più bassi, si possono notare coefficienti di autocorrelazione positivi rapidamente decrescenti e per i lag successivi  oscillazioni intorno allo zero. Ciò significa che esiste nella serie una persistenza di valori a breve termine, nel senso che se la grandezza in studio ha valore più elevato della media in un mese, lo sarà anche in uno o due mesi successivi e così per valori inferiori alla media.

-Se la serie storica presenta oscillazioni, anche il correlogramma tende ad assumere valori positivi e negativi, oscillando con lo stesso periodo della serie fino a smorzarsi ai lags più elevati. Se es. esiste un componente stagionale di periodo 12 mesi, nei dintorni del coefficiente di lag 12 ci sarà una zona significativamente diversa da zero.

– Nelle serie non stazionarie (presenza di TREND) i valori di rh non scendono velocemente a zero, ma si mantengono significativi per più valori del lag e solo se l’effetto del TREND è paragonabile alle altre eventuali relazioni presenti nei dati è possibile intuirle nel grafico (GRAFICO. N.2)

 

IL TEST DI DURBIN WATSON

Così la lettura dei correlogrammi talora può risultare ardua. Un modo veloce, affidabile e quantitativo per testare l’ipotesi che esista all’interno di una serie storica correlazione fra i suoi termini, cioè i termini non siano indipendenti, è somministrare alla serie il test di Durbin Watson ([8] 18-20), la cui statistica è espressa dalla formula:

 

d =∑ (ei – ei-1)2 /∑ ei2

 

La sommatoria al numeratore inizia dal 2° termine (i=2) e coinvolge ni termini . La statistica d varia da 0 a 4 e quando l’ipotesi nulla è vera (autocorrelazione assente) d dovrebbe essere vicino a 2. Il test permette di decidere di respingere l’ipotesi nulla, di accettarla o essere inconclusivo. Utilizzando la tabella opportuna   (allegata a queste note) si ottengono i valori critici di dl e du che servono per la decisione: all’interno dell’intervallo dl-du, la situazione è incerta; a sinistra di dl , si respinge l’ipotesi nulla. Vedremo in seguito come si calcola d con R e come si usa la tabella.

 

Il programma CORR, scritto in Qbasic, riportato in nota, permette il calcolo dei coefficienti di autocorrelazione con l’errore (un qualsiasi programma di grafica permetterà di costruire il correlogramma) e il calcolo della statistica di D. W.  Abbiamo già visto (vedere  programminosul correlogramma) come operare anche con il linguaggio R.

 

13 – APPENDICE2

APPENDICE2

Programmi in Qbasic e tabelle

PROGRAMMA CORR (coefficienti di autocorrelazione, il test di Durbin Watson, il test di Lin Mudholkar, Analisi spettrare per il periodogramma

 programma_period0001

 programma_period0002programma_period0003programma_period0004durbin watson_blogpag60001

i

durbin watson0002

i

durbin watson_blogpag60001

durbin watson0003

i

durbin watson0004

i

durbin watson0005

i

durbin watson0006

i

durbin watson0007

i

stat_reg_mlr_blog0001

durbin watson0008

GRAF. N.9

durbin watson0009

i

durbin watson0010

programma_period0005

i

durbin watson0003

durbin watson0001

14 – APPENDICE3

APPENDIX3

TABELLE  DEI RISULTATI

reg_tabelle_blog0001

i

reg_tabelle_blog0003

i

reg_tabelle_blog0004

i

reg_tabelle_blog0005

15 – APPENDICE4

APPENDIX4

ANALISI, CON IL LINGUAGGIO R, DELLA SERIE STORICA TRIMESTRALE RIVISITATA E AMPLIATA CON PERIODOGRAMMI E RISULTATI

period_det_trim2

FIG.1-20001

FIG.2'-3

FIG.4-6

 

FIG.8-11
FIG.7

 

FIG.120001

 

DA CAMBIARE:

> rm(list=ls(all=TRUE))
> #SCRIPTS IN R
>
> library(graphics)
> library(tseries)

‘tseries’ version: 0.10-32

‘tseries’ is a package for time series analysis and computational
finance.

See ‘library(help=”tseries”)’ for details.

> library(stats)
> #library(UsingR)
> library(lattice)
> library(lmtest)
Carico il pacchetto richiesto: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

as.Date, as.Date.numeric

>
> w=c(0.033,0.043,0.051,0.059,0.061,0.063,0.053,0.036,0.046,0.056,0.063,0.048,0.053,0.043,
+ 0.066,0.053,0.082,0.06,0.08,0.076,0.056,0.036,0.05,0.053,0.056,0.058,0.061,0.063,0.065,
+ 0.068,0.0815,0.095,0.079,0.063,0.069,0.074,
+ 0.08,0.0765,0.073,0.0695,0.066,0.093,0.083,
+ 0.073,0.063,0.074,0.067,0.06,0.086,0.08,0.073,0.067,0.089,0.064,0.087,0.079,0.07,0.065,0.06,.063)
>
> par(ask=T)
>
> par(mfrow=c(1,3))
>
> trim=matrix(w,ncol=3,byrow=T)
>
>
> medietrim=rowMeans(trim)
>
> medietrim
[1] 0.04233333 0.06100000 0.04500000 0.05566667 0.05400000 0.06500000
[7] 0.07066667 0.04633333 0.05833333 0.06533333 0.08516667 0.06866667
[13] 0.07650000 0.07616667 0.07300000 0.06700000 0.07966667 0.07333333
[19] 0.07866667 0.06266667
>
> # FIG.1
> ts.plot(medietrim,type=”l”,main=”FIG.1″) #finchè non lo sostituisco posso usare abline
Aspetto per confermare cambio pagina…
>
> w1=c(1:20)
> regtrim=lm(medietrim~w1)
> abline(regtrim)
>
> summary(regtrim)

Call:
lm(formula = medietrim ~ w1)

Residuals:
Min 1Q Median 3Q Max
-0.015979 -0.005078 0.001069 0.006031 0.019235

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0503921 0.0041790 12.058 4.67e-10 ***
w1 0.0014127 0.0003489 4.049 0.000752 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.008996 on 18 degrees of freedom
Multiple R-squared: 0.4767, Adjusted R-squared: 0.4476
F-statistic: 16.4 on 1 and 18 DF, p-value: 0.0007524

>
> val_pred_w=predict(regtrim) #calcolo i 20 valori predetti dalla prima regressione
> length(val_pred_w)
[1] 20
>
>
> detrend_trim=medietrim-val_pred_w
> detrend_trim
1 2 3 4 5
-0.0094714286 0.0077825815 -0.0096300752 -0.0003760652 -0.0034553885
6 7 8 9 10
0.0061319549 0.0103859649 -0.0153600251 -0.0047726817 0.0008146617
11 12 13 14 15
0.0192353383 0.0013226817 0.0077433584 0.0059973684 0.0014180451
16 17 18 19 20
-0.0059946115 0.0052593985 -0.0024865915 0.0014340852 -0.0159785714
>
> #FIG.2
> plot(detrend_trim,type=”l”, main=”FIG.2″)
>
>
> detrend_trim
1 2 3 4 5
-0.0094714286 0.0077825815 -0.0096300752 -0.0003760652 -0.0034553885
6 7 8 9 10
0.0061319549 0.0103859649 -0.0153600251 -0.0047726817 0.0008146617
11 12 13 14 15
0.0192353383 0.0013226817 0.0077433584 0.0059973684 0.0014180451
16 17 18 19 20
-0.0059946115 0.0052593985 -0.0024865915 0.0014340852 -0.0159785714
>
> trim1=matrix(detrend_trim,ncol=4,byrow=T)
> medietrim1=colMeans(trim1)
> medietrim1_5anni=rep(medietrim1,5)
>
> #FIG.3
> plot(medietrim1_5anni,type=”l”,main=”FIG.3″)
>
> medietrim1_5anni
[1] -0.0009393484 0.0036479950 0.0045686717 -0.0072773183 -0.0009393484
[6] 0.0036479950 0.0045686717 -0.0072773183 -0.0009393484 0.0036479950
[11] 0.0045686717 -0.0072773183 -0.0009393484 0.0036479950 0.0045686717
[16] -0.0072773183 -0.0009393484 0.0036479950 0.0045686717 -0.0072773183
>
> par(mfrow=c(2,2))
>
> #FIG.4
> acf(medietrim1_5anni,main=”FIG.4″)
Aspetto per confermare cambio pagina…
>
> valAdjtrim=medietrim-medietrim1_5anni #trend_ random
>
>
> fitadj_trim=lm(valAdjtrim~w1)
>
> fitadj_trim

Call:
lm(formula = valAdjtrim ~ w1)

Coefficients:
(Intercept) w1
0.049678 0.001481

>
> summary(fitadj_trim)

Call:
lm(formula = valAdjtrim ~ w1)

Residuals:
Min 1Q Median 3Q Max
-0.0136886 -0.0044597 -0.0006167 0.0058313 0.0146327

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.049678 0.003486 14.251 3.03e-11 ***
w1 0.001481 0.000291 5.088 7.67e-05 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.007504 on 18 degrees of freedom
Multiple R-squared: 0.5899, Adjusted R-squared: 0.5671
F-statistic: 25.89 on 1 and 18 DF, p-value: 7.671e-05

>
> #FIG.5
> plot(valAdjtrim,type=”l”,main=”FIG.5″)
> abline(fitadj_trim)
>
> #ANALISI RESIDUI
>
>
> dwtest(fitadj_trim, alternative=”two.sided”)

Durbin-Watson test

data: fitadj_trim
DW = 1.9024, p-value = 0.6301
alternative hypothesis: true autocorrelation is not 0

> #forse potremo interpolare l’elemento 11
>
> #FIG.6
> res=resid(fitadj_trim)
> plot(res,type=”l”, main=”FIG.6″)
>
> #FIG.7
> acf(res, main=”FIG.7″)
>
> par(mfrow=c(2,2))
> #FIG.8-12
> plot(fitadj_trim)
Aspetto per confermare cambio pagina…
>
>
>

Far girare il precedente programma. Applicare a detrend_trim il periodogramma  e trasformare in formula analitica l’oscillazione o le ascillazioni e provare a toglierla(toglierle) da medietrim (a da detrend _trim) per vedere se spariscono dal loro periodogramma i picchi rilevanti. E’ un buon metodo incrociato di testare il Periodogramma rivisitato.

16 – APPENDICE5

APPENDIX5

IL SENSO COMUNE, L’INSEGNAMENTO SCIENTIFICO ED I SAPERI PREPOSTI ALLE SCELTE – UN PRIMO APPROCCIO OPERATIVO ALL’ANALISI DI FOURIER COL SUPPORTO DEL COMPUTER  del dott. Piero Pistoia

0 – LA PREMESSA

MATH_FOURIER_PREMESSA1

BIBLIOGRAFIA DELLA PREMESSA

four_bibl0001

 

1 – L’ARTICOLO GUIDA

L’articolo sull’analisi di fourier su dati reali e simulati col Mathematica di Wolfram 4.2                 dott. Piero Pistoia

ARTFOUART-math

Seguirà la diretta trascrizione

 2 – IL PROGRAMMA CON ESERCITAZIONI

Analisi di serie storiche reali e simulate dott. Piero Pistoia

ATTENZIONE! le linee di programma attive non sono incluse fra apici. Cambiando opportunamente le inclusioni di linee nei diversi segmenti del programma, si possono fa girare i diversi esempi, e proporne di nuovi.

A0-Esempio N.0

ANALISI DELL’ESEMPIO N° 0

Le linee attive di questo esempio sono state evidenziate

Si trascriva manualmente o con copia/incolla i seguenti scripts  sulla consolle del MATHEMATICA DI WOLFRAM (vers. non superiore alla 6.0); si evidenzi e si batta shift-enter: si otterranno i risultati ed i grafici non inseriti.

"Si forniscono diversi vettori di dati sperimentali di esempio immessi 
    direttamente o tramite Table; per renderli attivi basta eliminare agli 
    estremi le virgolette.Se l'analisi diventa più complessa rispetto ad una 
ricerca di armoniche di Fourier a confronto con la serie iniziale, si può 
utilizzare il SEGMENTO DELLE REGRESSIONI (lineare e quadratica) per ottenere 
yg2 ed il SEGMENTO DELLE ARMONICHE RILEVANTI (yg3 e yg4) individuate in una 
prova precedente. Abbiamo da sostituire il nome di qualche vettore e aprire o 
chiudere (cancellando o inserendo virgolette) istruzioni nei diversi segmenti 
del programma secondo ciò che vogliamo fare. In yg1 c'è il vettore dati 
iniziale. In yg2 c'è il vettore detrendizzato. In yg3, quello delle armoniche 
rilevanti. In v, il vettore di Fourier fornito dall'analisi. Altri segmenti 
su cui intervenire: IL GRAFICO ygf dove va inserita la variabile (ygi) da 
confrontare con Fourier (v); il segmento di IMPOSIZIONE NUMERO ARMONICHE m; 
il segmento di SCELTA VARIABILI DA SOTTOPORRE A FOURIER (ygi); il segmento 
per cambiare la variabile nell'ERRORE STANDARD. In ogni esempio si accenna 
alle modifiche specifiche da apportare ai diversi segmenti";


"ESEMPIO N.0";
"Esempio illustrativo riportato alle pagine 3-4 dell'art.: imporre il numero \
di armoniche m=1 oppure 2 nel segmento relativo e confrontare il grafico ygf \
che gestisce la variabile yg1 dei dati seguenti, con quello di v (ygf1); \
controllare infine i risultati con i dati del testo";

yt=N[{103.585, 99.768, 97.968, 99.725, 101.379, 99.595, 96.376, 96.469, \
100.693, 104.443}]

"ESEMPIO N.1"
"Si sottopongono a Fourier i dati tabellati seguenti(yg1). Si confrontano yg1 \
(tramite ygf) e v di Fourier(tramite ygf1); calcolo automatico di m";
 "yt=N[Table[100+4 Sin[2 t/21 2 Pi-Pi/2]+3 Sin[4 t/21 2 Pi+0]+6 Sin[5 t/21 2 \
Pi-1.745], {t,21}]]"

"Si detrendizza yg1 seguente ottenendo yg2 (si liberino le istruzioni del \
segmento TREND), che poniamo come variabile in ygf (si inserisca nella sua \
espressione); si sottopone yg2 a Fourier (v) nel segmento "SCELTA VETTORE \
DATI"; confrontiamo ygf1 (grafico di v) e ygf; inserire variabile yg2 \
nell'espressione errore"  
"yt=N[Table[100+4 Sin[2 t/21 2 Pi-Pi/2]+3 Sin[4 t/21 2 Pi+0]+6 Sin[5 t/21 2 \
Pi-1.745]+0.5 t + (Random[]-1/2),{t,21}]]";

"ESEMPIO N.2";
"Si utilizza il vettore originale yg1 e si confronta con v di Fourier (ygf \
con ygf1), come nell'esempio N.1, prima parte; m automatico. Esempio \
interessante per controllare come Fourier legge i dati"
"yt=Table[N[Sin[2 Pi 30 t/256]+.05t+(Random[]-1/2)],{t,256}]"

"ESEMPIO N.3";
"Come l'esempio N.2. Ci insegna come Fourier <sente> i coseni"
"yt=N[Table[100+4 Cos[2 t/21 2 Pi-Pi/2]+3 Cos[4 t/21 2 Pi+0]+6 Cos[5 t/21 2 \
Pi-1.745],{t,21}]]";

"ESEMPIO N.4";
"Come il N.2. Ci assicura del funzionamento del programma"
"yt= N[Table[100+3 Sin[2 Pi 2 t/21+ Pi/2],{t,1,21}]]";
"I dati successivi sono stati campionati da Makridakis combinando l'esempio \
precedente ed un random (pag. 402, [])";
"yt={106.578,92.597,99.899,97.132,93.121,95.081,102.807,106.944,100.443,95.\
546,103.836,107.576,104.658,91.562,91.661,97.984,111.195,100.127,94.815,105.\
009,110.425}";

"ESEMPIO N.5"
"Prima parte.
Si detrendizza il vettore dati yg1, liberando, nel segmento TREND, il calcolo \
dei coefficienti B0 e B1 della retta interpolante, trovando yg2 che \
inseriremo in ygf nel segmento GRAFICO DA CONFRONTARE CON FOURIER. Calcolo \
automatico di m. Nel segmento SCELTA VETTORE PER FOURIER, poniamo yg2 in yt e \
nella formula dell'ERRORE STANDARD. Si fa girare il programma una prima volta \
e si osservano quali sono le armoniche rilevanti. Di esse si ricopiano i  \
parametri trovati (numero armonica, ampiezza, fase), con i quali  si \
tabellano le 4 armoniche rilevanti, trascrivendole nel segmento ARMONICHE \
RILEVANTI.
Seconda parte.
Nel segmento ARMONICHE RILEVANTI si tabellano le espressioni di queste 4 \
armoniche, si sommano i relativi vettori in yg3. Si liberano queste 4 \
armoniche e il loro vettore somma yg3. Si pone poi la variabile yg3 in ygf \
per confrontare yg3 con v (risultato di Fourier su yg2). Calcolo automatico \
di m. Si sceglie per Fourier la variabile yt=yg2. Nel segmento dell'ERRORE si \
pone yg3 e si rilancia il programma una seconda volta. E' un modo per \
cogliere le uniformità periodiche all'interno di dati storici"

"yt=N[{0.0330,0.0430,.0510,.0590,.0610,.0630,.053,.036,.0460,.0560,.0630,.\
0480,.0530,.0430,.0660,.053,.0820,.0600,.0800,.0760,.0560,
.0360,.0500,.053,.0560,.0580,.0610,.0630,.0650,.0680,.0815,.095,
.0790,.0630,.0690,.0740,.0800,.0765,.0730,.0695,.0660,.0930,.0830,
.0730,.0630,.0740,.0670,.06,.0860,.0800,.0730,.0670,.0890,.0640,
.0870,.079,.0700,.0650,.0600,.0630}]"

"Le successive righe sempre attive"
yg1=yt
n=Length[yt];

"SEGMENTO DELLE REGRESSIONI"

"f[x_]:=Fit[yt,{1,x,x^2},x]"
"f[x_]:=Fit[yt,{1,x},x]"

"yt1=N[Table[f[t],{t,60}]]?"
"La precedente istruzione dà problemi"

"Trovo l'ordinata all'origine e la pendenza"
"B0=f[x]/.x\[Rule]0"
"f1=f[x]/.x\[Rule]1"
"B1=f1-B0"
"B2=B0+B1 t"

"Un secondo modo di trovare B0 e B1";

"xt=Table[i, {i, 1, n}]";
"a=xt yt";
"Sxy=Apply[Plus, xt yt]";
"Sx=Apply[Plus, xt]";
"Sy=Apply[Plus, yt]";
"xq=xt^2";
"Sxq=Apply[Plus, xq]";
"yq=yt^2";
"qSx=Sx^2";
"B1=(n Sxy-Sx Sy)/(n Sxq-qSx)";
"B0=Sy/n-B1 Sx/n";
"B2=B0+B1 t";

"Tabello la retta"

"yt1=N[Table[B2,{t,n}]]"
"yt1=Flatten[yt1]"
"In yt1 ci sono i dati relativi alla retta di regressione"
"In yg2 c'è il vettore detrendizzato dei dati iniziali"
"yg2=yt-yt1"

"SEGMENTO DELLE ARMONICHE RILEVANTI"

"y4=Table[N[.004(Sin[.1333 Pi t+6.266])],{t,60}]";
"y5=Table[N[.007(Sin[.1667 Pi t+4.782])],{t,60}]";
"y8=Table[N[.004(Sin[.2667 Pi t+4.712])],{t,60}]";
"y9=Table[N[.004(Sin[.3000 Pi t+3.770])],{t,60}]";

"yg3=N[y4+y5+y8+y9]";

"In yg3 c'è il vettore dati di tutte le armoniche considerate rilevanti da \
precedente analisi. Se il programma passa da questo punto,
    ha senso misurare per es. la differenza con il vettore di tutte le \
armoniche di Fourier sui dati detrendizzati yg2";

"yg4= N[yg3+yt1]";

"In yg4 c'è il vettore di tutte le componenti considerate rilevanti compreso \
il trend. Ha senso un confronto fra i dati iniziali Yg1 o v (vettore di \
Fourier) e Yg4 "

" IL GRAFICO ygf E' DA CONFRONTARE CON QUELLO DI FOURIER ygf1"
" La variabile nel ListPlot successivo rappresenta il vettore da confrontare \
con la combinazione di armoniche di Fourier applicato ad un vettore dati. yg \
rappresenta il grafico di tale vettore"
"In ygi (i=1,2...) ci va il vettore da confrontare con v"

ygf=ListPlot[yg1,PlotJoined\[Rule]True,PlotRange\[Rule]Automatic,
    				   GridLines\[Rule]{Automatic,Automatic},
    AxesLabel\[Rule]{"Tempo","Dati \
(unità)"},PlotLabel\[Rule]FontForm["DOMINIO  DEL TEMPO",{"Times",12}]]

"CALCOLO AUTOMATICO DEL NUMERO ARMONICHE"
ny=Length[yg1]
n=ny;m=Mod[n,2]
If[m>0,  m=(n-1)/2, m=n/2-1]
"IMPOSIZIONE MANUALE NUMERO ARMONICHE"
m=1
"m=2"

"SCELTA VETTORE DATI DA SOTTOPORRE A FOURIER"
"IN yt CI SONO I DATI CHE VOGLIO ANALIZZARE CON FOURIER E L'ANALISI E' POSTA \
IN v"
"Se voglio analizzare con Fourier i dati iniziali:"
yt=yg1
"Se voglio analizzare i dati detrendizzati:"
"yt=yg2"
"Se voglio analizzare i dati relativi alle armoniche considerate rilevanti:"
"yt=yg3"
"Se voglio analizzare i dati di tutte le componenti rilevanti:"
"yt=yg4"

"VALORI DEL PARAMETRO ak="

"Calcolo gli ak con il comando Sum, sommando cioè gli n prodotti yt * la \
funzione coseno, per t=1 a n; faccio questo per ogni valore di k (da k=0 a \
n/2)tramite Table"

a1=Table [Sum[yt[[t]] Cos[2 Pi k t/n],{t,1,n}],{k,0,m}];
a=2*a1/n;

"Divido per due il primo elemento, per ottenere ao=media; Sopprimo poi il \
primo elemento"
a0=a[[1]]/2
a=Delete[a,1]
a=Chop[%]

"VALORI DEL PARAMETRO bk="

"Calcolo ora bk con la funzione seno con lo stesso procedimento di ak"
b1=Table[Sum[yt[[t]] Sin[2 Pi i t/n],{t,1,n}],{i,1,m}];
b=2 b1/n
b=Chop[%]
"Mentre ao/2 rappresenta la media, bo è sempre nullo"
b0=0

"AMPIEZZE ="

"Con ak e bk calcolo le ampiezze e le fasi dell'f(t) iniziale; Individuo il 
vettore dei numeri da mettere sulle ascisse nel dominio della frequenza 
(i/n o n/i) e con i vettori xi e yi 
costruisco la lista {xi,yi}; disegno infine i plots"
ro=Sqrt[a^2+b^2]
ro=N[Chop[%]]
ro=Flatten[ro]


Theta={}
i=1
While[i<m+1,
    f2=Abs[a[[i]]/b[[i]]];
    f2=180/Pi ArcTan[f2];
    If[b[[i]]>0 && a[[i]]>0 , Theta=N[Append[Theta,f2]]];
    If[b[[i]]<0 && a[[i]]>0, Theta=N[Append[Theta,180-f2]]];
    If[b[[i]]<0 && a[[i]]<0, Theta=N[Append[Theta,180+f2]]];If[b[[i]]>0 && a[[
    i]]<0, Theta=N[Append[Theta,360-f2]]];
    If [(a[[i]]==0 && b[[i]]==0),Theta=N[Append[Theta,0]]]; 
     If[((
    b[[i]]<0 || b[[i]]>0) && a[[i]]\[Equal]0),Theta=N[Append[Theta,0]]];
    
     If[b[[i]]\[Equal]0 && a[[i]]>0 ,Theta=N[Append[Theta,90]]];
    If[b[[i]]\[Equal]0 && a[[i]]<0, Theta=N[Append[Theta,-90+360]]]; i++];

"FASE ="

Theta=Theta

"Theta=N[ArcTan[a,b]*180/Pi]"

"RISULTATI DI FOURIER"
v=Table[a0+Sum[(a[[k]] Cos[2 Pi k t/n]+b[[k]] Sin[2Pi  k \
t/n]),{k,1,m}],{t,1,n}];

"GRAFICO RISULTATI DI FOURIER (ygf1)"
ygf1=ListPlot[v,PlotJoined\[Rule]False,GridLines\[Rule]{Automatic,Automatic},
    PlotLabel\[Rule]FontForm["GRAFICI  DI CONTROLLO",{"Times",12}]]

"CONFRONTO"
ygf2=Show[ygf,ygf1,PlotRange\[Rule]{Automatic,Automatic}]

"Calcolo l'ERRORE STANDARD DELLA STIMA"

ESS=Sqrt[Apply[Plus,(yg1-v)^2]/(n-2)]

"x=N[Table[i,{i,1,m}]]";
"c1=x";
Length[x];
Length[ro];
"For[i=1,i<m,i++,j=i*2;c=c1;yi=ro[[i]];
  c=Insert[c,yi,j];c1=c]";
"d1=Partition[c,2]";
Needs["Graphics`Graphics`"]
BarChart[ro]
ListPlot[ro, PlotJoined\[Rule]True,PlotRange\[Rule]All, 
  GridLines\[Rule]{Automatic,Automatic},AxesOrigin\[Rule]{0,
          0},AxesLabel\[Rule]{"Cicli in n dati", "Ampiezza \
"},PlotLabel\[Rule]FontForm["DOMINIO  DELLA FREQUENZA",{"Times",12}]]
"c1=x";
"For[i=1, i<m,i++,j=i*2;c=c1;yi=Theta[[i]];
  c=Insert[c,yi,j];c1=c]";
"d2=Partition[c,2]";
ListPlot[Theta, PlotJoined\[Rule]True,PlotRange\[Rule]All, \
GridLines\[Rule]{Automatic,Automatic},AxesOrigin\[Rule]{0,0},
  AxesLabel\[Rule]{"Frequenza","Fase"}]

———————————————————————————–

A1-Esempio N.1

————————————————————————

A2-Esempio N.2

math_es_20001

math_es_20002

math_es_20003



math_es_20004

math_es_20005

math_es_20006

math_es_20007RISULTATI ESEMPIO 2

math_es_20008
math_es_20009

math_es_20010

math_es_20011

math_es_130001math_es_20011
math_es_130001
math_es_150001

math_es_20016

A5-Esempio N.5

 

Serie detrendizzata delle concentrazioni As 

ANALISI DEI DATI REALI DELL’ESEMPIO N° 5

priodogramma0001

L’IDEA E’ QUESTA:

– SUI SESSANTA DATI DELLA CONCENTRAZIONE ARSENICO (yt, GRAF. N.1) IN ALCUNE SORGENTI DELLA CARLINA (PROV. SIENA), SI FA UNA REGRESSIONE LINEARE ED I SUOI  60 VALORI PREDETTI  SI SOTTRAGGONO DA yt, OTTENENDO LA SERIE DETRENDIZZATA.

– QUEST’ULTIMA SI SOTTOPONE AL PERIODOGRAMMA CHE, IN USCITA, PERMETTE DI CALCOLARE LE SUE COMPONENTI ARMONICHE.

– SOMMANDO LE COMPONENTI ARMONICHE RILEVANTI PIU’ I VALORI DEL TREND E SOTTRAENDO TALE SOMMA DALLA SERIE ORIGINALE yt, SI OTTERRA’ LA “STIMA DELL’ERRORE STANDARD” CHE DA’ UN’IDEA DELLA BONTA’ DEL PROCESSO.

Si trascriva manualmente i seguenti scripts  sulla consolle del MATHEMATICA DI WOLFRAM (vers. non superiore alla 6.0); si evidenzi e si batta shift-enter: si otterranno i risultati e grafici riportati alla fine di questo programma (si noti in particolare il grafico ampiezza-numero armoniche eseguito sulla serie detrendizzata, dove è evidente il picco all’armonica n°5)

Period_con_math0001

Period_con_math0002

period_con_math10001

period_con_math10002

period_con_math10003

——————————————————————–

RISULTATI DEL PROGRAMMA ESEMPIO N.5 (conc. As detrend)

Period_con_math20001

Period_con_math20002

Period_con_math20003x

Period_con_math20004

Period_con_math20005x

Period_con_math20006

Period_con_math20007x

A4-Esempio N.4

ANALISI DELL’ESEMPIO N° 4 CON RISULTATI E GRAFICI (DATI SIMULATI)
Mathematica0001

Mathematica0002

Mathematica0003

Mathematica0004 - Copia

Mathematica0005 - Copia

Mathematica0006

Mathematica10001

Mathematica10002

Mathematica10003

Mathematica10004

Mathematica20001 - Copia

Mathematica20002 - Copia

Mathematica20003

Mathematica20004

A6-Esempio N.6

Oscillazione mensile ozono a Montecerboli (Pomarance,Pi); 2007-2011

fouroz20001

fouroz20002

fouroz20004

fouroz20006

fouroz20005

fouroz20006

fouroz20007fouroz2008RISULTATI GRAFICI OZONO
fouroz2009

fouroz20010

———————————————————————————-

 

ESEMPIO N° 5 CHE USA LE ARMONICHE RILEVANTI MESSE IN FORMULA IN UNA PRE-PROVA

four_art1 da correggere

four_art2

four_art3

four_art4

four_art5

four_art6

four_art7

four_art8

four_art9

four_art10

four_art11

four_art12

four_art13

four_art14

four_art15

four_art16

four_art17

four_art18

fuor_art19 da correggere

Annunci

UN PERCORSO VERSO IL ‘PERIODOGRAMMA’ ED ALTRO: CONCETTI PORTANTI DELL’ANALISI DELLE SERIE STORICHE a cura di Pf. Bianchi – P. Pistoia

Curriculum di piero pistoia:

PIERO PISTOIA CURRICULUMOK

 

CHI E’ L’AUTORE (traccia): 

Piero Pistoia, diplomato negli anni ’50 presso il Liceo Classico Galileo Galilei di Pisa, è dottore in Scienze Geologiche con lode e, da borsista, ha lavorato e pubblicato presso l’Istituto di Geologia Nucleare di Pisa, misurando le età degli “strani” graniti associati alle ofioliti (1) e studiando i serbatoi di gas e vapori della zona di Larderello. Successivamente ha scritto una cinquantina di articoli pubblicati a stampa, a taglio didattico-epistemologico, di cui circa la metà retribuiti secondo legge,  dagli editori Loescher, Torino, (rivista “La Ricerca”), La Scuola di Brescia (“Didattica delle Scienze”), a controllo accademico ed altri, affrontando svariati problemi su temi scientifici: dall’astrofisica all’informatica, dall’antropologia culturale all’evoluzione dell’uomo, dalla fisica alla matematica applicata e alla statistica, dalla geologia applicata al Neoautoctono toscano, dall’origine dell’Appennino alla storia delle ofioliti, alle mineralizzazioni delle antiche cave in Val di Cecina (in particolare su calcedonio, opale e magnesite) ecc..  En passant, ha scritto qualcosa anche sul rapporto Scienza e Poesia, sul perché la Poesia ‘vera’ ha vita infinita (per mere ragioni logiche o perché coglie l’archetipo evolutivo profondo dell’umanità?); ha scritto alcuni commenti a poesie riprese da antologie scolastiche e,  infine decine di ‘tentativi’ poetici senza pretese. Molti di tali lavori sono stati riportati su questo blog. (2)
NOTE
(1) L’età dei graniti delle Argille Scagliose, associati alle ofioliti, corroborò sia l’ipotesi che esse fossero ‘strappate’ dal basamento ercinico durante i complessi  eventi che costruirono la catena appenninica, sia, indirettamente, rafforzò la teoria a falde si ricoprimento nell’orogenesi appenninica. Fu escluso così che il granito associato alle ofioliti derivasse, almeno non in tutti i casi, da una cristallizzazione frazionata (serie di Bowen) da un magma basico od ultrabasico.
(2) Piero Pistoia ha superato concorsi abilitativi nazionali, al tempo fortemente selettivi (cioè non frequentò mai i famigerati Corsi Abilitanti, fortemente voluti dai sindacati dei docenti!), per l’insegnamento in particolare nella Scuola Superiore per le seguenti discipline: Scienze Naturali, Chimica, Geografia, Merceologia, Agraria, FISICA e MATEMATICA. Le due ultime materie sono maiuscole per indicare che Piero Pistoia in esse, in tempi diversi, fu nominato in ruolo, scegliendo poi la FISICA, che insegnò praticamente per tutta la sua vita operativa.

Pochi anni prima che l’ITIS di Pomarance fosse aggregato al Commerciale di Volterra, il dott. prof. Piero Pistoia fu nominato Preside Incaricato dal Provveditorato agli studi di Pisa, ottenendo il massimo dei voti.

 

 

ATTENZIONE QUESTO LAVORO VIENE CONTINUAMENTE DANNEGGIATO SPECIALMENTE NEGLI SCRIPTS NONOSTANTE LE CORREZIONI! (Piero Pistoia)

UN PERCORSO VERSO IL ‘PERIODOGRAMMA’ 
ATTRAVERSO ALCUNI CONCETTI PORTANTI DELL’ANALISI DELLE SERIE STORICHE

ATTENZIONE: INTERVENTI IN VIA DI COSTRUZIONE! La via si fa con l’andare!

Noi autori. di questo scritto, prof. Piero Pistoia, docente di Fisica, ed il collega di Matematica, prof. Pier Francesco Bianchi, stiamo lavorando ad un progetto per cercare di automatizzare il processo di analisi di una serie storica reale, usando il linguaggio R ed il MATHEMATICA di Wolfram, cercando di usare comandi di basso livello. Ci apriamo il sentiero nell’andare per cui non sarà rettilineo ma complesso e intrecciato e spesso dovremo tornare indietro, persi nei meandri di questo ‘zibaldone’  alla Leopardi di statistiche e programmazioni. Una cosa è certa però: ci stiamo divertendo e andremo avanti e forse chissà se ci fermeremo.

L’obbiettivo di questo lavoro è fornire strumenti direttamente operativi al lettore perchè possa scoprire all’interno di una serie storica componenti rilevanti della sua variazione. Faremo una sintesi argomentativa sui concetti di statistica implicati nello studio di una serie storia (correlogramma, periodogramma, media mobile, regressione lineare, analisi dei residui…) fino a ‘leggere’ all’interno di serie specifiche, seguendo un itinerario guidato.

Abbiamo iniziato “forzando” R a raccogliere informazioni preliminari sulla serie in studio per poi intervenire su essa correggendo ed aggiustando ( correzione degli outliers, interpolazioni per mancanza di dati, ecc.). Sorge poi il problema di cosa cercare all’interno dei dati e come cercarlo a partire da grafici e altri tests (correlogrammi e periodogrammi ecc.). A questo punto si è posto il problema di “fittare” una combinazione di onde del seno ai dati di una serie storica sperimentale usando il metodo dei minimi quadrati. Una volta imparato come calcolare i coefficienti di una regressione lineare, appunto condotta sui dati con una combinazione polinomiale di funzioni sinusoidali, cercheremo di precisare con esempi il processo chiamato Analisi armonica di Fourier.

Ma oggi, in questo contesto, ci siamo soffermati a riflettere su come possiamo, utilizzando sempre il linguaggio R e il MATHEMATICA di Wolfram, costruire il cosiddetto Periodogramma di una serie storica che possa individuare, se ci sono, oscillazioni sinusoidali rilevanti. Il Periodogramma in generale riporta sulle ascisse valori in successione, ognuno dei quali rappresenta il numero delle oscillazioni complete che quella particolare onda compie in n dati (così,se n=60, 5 rappresenta la presenza nei dati di una oscillazione che fa 5 cicli completi in 60 dati, cioè di periodo 12). Le oscillazioni rilevanti, suggerite dal Periodogramma, possono essere tolte poi dai dati originali (destagionalizzazione), per esempio, con medie mobili opportune, onde iniziare la scomposizione della serie stessa e così via.
Vi proponiamo allora, nel proseguio, attraverso un percorso piuttosto lungo ed articolato, punteggiato di concetti ed esempi di calcolo, come obbiettivo ‘regolativo’ sl, nostre versioni di Periodogramma, da noi scritte in R e in Mathematica di Wolfram, che pur non ottimizzate, riescano efficacemente a funzionare su due esempi, uno da noi controllato ed l’altro reale.

Tutto quello che verrà detto durante questo lavoro, pur non avendo la pretesa di esaurire le problematiche ivi implicate (per questo vedere bibliografia), speriamo aiuterà il lettore, se interessato, a seguire l’analisi di una serie storica, attraverso cammini meno usuali e teorici, fornendo strumenti operativi per poter affrontare studi più organizzati in un secondo momento.

I conti possono essere seguiti su una qualsiasi spread-sheet oppure attraverso tre programmi  in Qbasic allegati (scritti da P. Pistoia), poco curati nella forma, ma che contengono routines efficaci, e/o utilizzando, come abbiamo accennato, i comandi di due grossi programmi di statistica, il programma R ed il linguaggio del Mathematica di Wolfram.

IL CORRELOGRAMMA: UNA ‘LENTE’ PER OSSERVARE IN UNA SERIE STORICA VARIE COMPONENTI DELLA SUA VARIAZIONE

All’interno di una serie storica reale, per es., dei 60 ( mg/l) dati mensili yt della concentrazione As nelle acque della Carlina, che a più riprese analizzeremo, è possibile individuare varie serie componenti: a) una o più serie stagionali, nel senso lato della parola, cioè oscillazioni regolari, periodiche, in accordo col calendario o con l’orologio; b) una serie corrispondente al trend o variazione a lungo termine della media, che di fatto comprende tutte le componenti cicliche la cui lunghezza d’onda supera la serie temporale osservata, cioè supera lo stesso range dei dati; c) uno o più cicli, cioè fluttuazioni irregolari intorno ad una tendenza a lungo termine, non periodiche (almeno all’interno dei dati) talora intermittenti, di durata variabile, non correlate al calendario, i cui effetti durano oltre i dati, ma coglibili anche all’interno del range dei dati,  d) la componente irregolare o casuale (random) che riassume lo “white noise” e infine e) qualche componente occasionale o erratica (periodo di siccità, terremoto, sciopero…).

Sintetizzando, possiamo ragionevolmente affermare che i dati annuali tipicamente esibiscono un trend ed un effetto ciclico, ma non un effetto stagionale. Ci sono però anche eccezioni: si pensi agli effetti stagionali di periodo superiore a qualche anno di alcuni fenomeni astronomici (ciclo delle macchie solari, di periodo circa 11 anni, che forse può attivare altri fenomeni analoghi nell’atmosfera). I dati mensili  e quelli ‘quaterly’ (relativi ad 1/4  di anno, cioè trimestrali), probabilmente mostrano trend e influenze stagionali, ma usualmente non componenti cicliche a meno che i dati coinvolgano molti anni. In generale gli effetti stagionali sono visti come oscillazioni ripetitive entro il tempo di un anno. Sono comprese negli effetti stagionali anche le oscillazioni all’interno delle 24 ore per dati orari e le oscillazioni all’interno della settimana e all’interno del mese per dati giornalieri.

Prima di cercare di scomporre nelle sue componenti elementari una serie storica, conviene dare uno sguardo al suo contenut con un grafico opportuno,  con lo strumento Correlogramma ( [1] pagg. 376-390)  e  con il Periodogramma, che presenteremo definitivo al termine del lavoro, strumenti essenziali anche per l’analisi dei RESIDUI.

Il coefficiente di correlazione di Pearson misura la correlazione fra due variabili aleatorie che dipendono linearmente l’una dall’altra. Ha valore +1 se le due variabili variano linearmente in fase e -1 se variano in controvase. Questo coefficiente acquista valori fra -1 e +1. Il coefficiente di correlazione di Pearson dipende dalla covarianza, covarianza=Somma[(x-xm)*(y-ym]/n, e dalla deviazione standard di entrambe le variabili: coeff. di correlazione=Covarianza/(DSx*DSy).

I coefficienti di auto-correlazione rh, dove h=0,1,2…q con q minore od uguale a (n-2)/2, sono coefficienti di correlazione, calcolati per ogni valore di h, che misurano la concordanza o la discordanza fra i valori di una serie storica e quelli della stessa però slittati di h unità di tempo (lag h), consentendo di analizzare la sua struttura interna, ossia i legami fra i termini della stessa ([8] 18-20).

rh = Σi[(y(t)-ym)(y(t+h)-ym)]/[(n-h)*Σj(y(t)-ym)^2/n)] dove i va da t=1…n-h e j va da t=1 … n-h # da togliere l’ultima h!!!

in alcuni testi viene abolito il fattore n/(n-h).

Tale formula presenta la semplificazione di poter utilizzare una media unica per le yt (quella dei dati originali), presupponendo una situazione stazionaria ([8] pag. 19 e [2] pag. 133). In particolare ro=1 (lag h=0), nessun slittamento e gli altri rh assumono valori fra +1 (completa concordanza) e -1 (totale discordanza).
Il correlogramma è la rappresentazione grafica dei coefficienti di auto-correlazione, che sono (n-2)/2, in funzione degli slittamenti (lag h) e ci permette di vedere se la serie storica possiede qualche regolarità interna.
I coefficienti di auto-correlazione di dati random hanno distribuzione campionaria che può essere approssimata da una curva gaussiana con media zero e errore standard 1/√(n). Questo significa che il 95% di tutti i coeff. di auto-cor., calcolati da tutti i possibili campioni estratti, dovrebbero giacere entro un range specificato da: zero ± 1.96*1/√(n) (1.96 errori standard). I dati cioè della serie saranno da considerare random, se i coeff. di auto-cor. staranno entro i limiti:

-1.96*(1/√n) ≤ rh ≤ +1.96*(1/√n)      fascia dell’errore: +/- 2/√n

Per l’interpretazione dei correlogrammi vedere [8] 20-25.

– In una serie storica completamente casuale, i cui successivi valori sono considerati tutti indipendenti fra loro (non correlati), tutti i valori rh,  eccetto ro che è sempre +1 ( correlazione della serie con se stessa), oscilleranno  in maniera casuale intorno allo zero entro la fascia dell’errore.

– I coefficienti di correlazione per dati stazionari (assenza di trend) vanno velocemente a zero dopo 3-4  lags di tempo (forse fino a 5), mentre nella serie non stazionaria, essi sono significativamente diversi da zero per varie unita di tempo anche se tendono a diminuire, es., serie che contiene trend (vedere graf. yt). Nella serie stazionaria esiste una persistenza di valori positivi o negativi a breve termine (se per es., il valore è più alto della media in un mese, lo sarà anche in uno o due mesi successivi. Data la brevità di questo fenomeno (fino a 5 lags max) si riscontra anche in correlogrammi di componenti erratiche.

-Anche col periodogramma (talora detto spettrogramma) è possibile individuare componenti oscillanti, ma anche i trends così da poterli eliminare dalla serie. Come primo input: qualsiasi serie temporale composta di n osservazioni ugualmente spaziate può essere decomposta tramite il metodo dei minimi quadrati in un numero di onde del seno di data frequenza, ampiezza e fase, soggette alle seguenti restrizioni: se n è dispari, allora il numero max di onde fittate è (n-1)/2, se n è pari tale numero è N/2-1. Da notare che in una serie temporale discreta, poichè non appaiono di fatto angoli da trattare, nella definizione di lunghezza d’onda e fase si fa riferimento alle unità di tempo usate per definire la serie, o al numero delle osservazioni n che ‘costruiscono’ la lunghezza d’onda.

– Se la serie storica presenta oscillazioni (stagionali:  oscillazioni orarie, giornaliere, settimanali, mensili …), anche il correlogramma tende ad assumere valori positivi e negativi, oscillando con lo stesso periodo della serie in studio fino a smorzarsi ai lags più elevati. Inoltre, se esiste, per es., una componente stagionale di periodo 12 mesi, il valore corrispondente al lag 12 sarà significativamente diverso da zero. 

Talora però la lettura del correlogramma risulta ardua. Un modo veloce e quantitativo per testare l’ipotesi che esista all’interno di una serie storica correlazione fra i suoi termini, cioè i termini non siano indipendenti, è somministrare alla serie il test di Durbin e Watson ([3] 949-951), la cui statistica è espressa dalla formula:

d = Σ[e(i)-e(i-1)^2]/Σei^2

La sommatoria al numeratore inizia al 2° termine (i=2= e coinvolge n-1 termini. la statistica d varia da 0 a 4 e quando l’ipotesi nulla è vera (autocorrelazione assente) d dovrebbe essere vicino a 2. Il test permette di decidere di respingere l’ipotesi nulla, di accettarla ovvero essere inconclusivo. Utilizzando la tabella opportuna  (allegata) si ottengono i valori critici dl e du che servono per la decisione: all’interno dell’intervallo dl-du, la situazione è incerta; a sinistra di dl, si respinge  l’ipotesi nulla.

Il programma CORR, scritto da P. Pistoia nel glorioso Qbasic (saturo di nostalgia e giovinezza) e allegato a questo scritto, permette il calcolo dei coefficienti di autocorrelazione con l’errore e il calcolo della statistica di Durbin-Wtson; un qualsiasi programma di grafica  poi permetterà di costruire il correlogramma sottoponendo ad esso i coefficienti di auto-corr. trovati. Il correlogramma verrà naturalmente ottenuto anche direttamente con i linguaggi di livello R ed il Mathematica di Wolfram. CORR infine fa anche l’analisi armonica.

ESERCIZI DI COSTRUZIONE E ‘LETTURA’ DEL CORRELOGRAMMA

Consideriamo la serie storica mensile yt di 60 dati della concentrazione As già nominata. Inseriamo in CORR vari vettori dati da analizzare, mutuati da questo studio (per es., yt, ESAs (Effetto stagionale), Yt1=yt-ESAs (Ciclo+Trend+Random), Yt1_smussato dai random …), nei comandi DATA al capoverso SERIE STORICA ORIGINALE. Successivamente possiamo trasformare in remarks le precedenti linee yt e aprire quelle successive, EFFETTO STAGIONALE As (ESAs), se gira il programma, troveremo il dati del correlogramma per la serie ESAs. Per aprire e chiudere i DATA si usa l’apice . Naturalmente possiamo introdurre la serie che vogliamo (es., RESIDUI..). Il programma nei due casi accennati, lanciato, fornirà 1) La tabella dei coefficienti di auto-correlazione 2)La statistica di Durbin Watson per controllare se c’è autocorrelazione nella serie in studio 3)La statistica di LinMudholcar che è un test sulla gaussiana; 4) L’analisi di Fourier (per ora sospesa). Basta riscrivere le linee di programma nella console del Qbasic.
Scriveremo poi un programmino in R per ottenere gli stessi risultati, molto più agile per avere comandi di livello superiore .Chi vuole può ottenere gli stessi risultati anche con lo straordinario Mathematica di Wolfram o addirittura con EXCEL.

SCRIPT IN Qbasic

‘ PROGRAMMA N. 1

CLS
PRINT TAB(15); “TESTS DI AUTOCORRELAZIONE E NORMALITA'”
PRINT
PRINT TAB(14); “Programma scritto a cura del DOTT. PIERO PISTOIA”
LOCATE 10, 5
PRINT “Il programma calcola: “: PRINT
PRINT TAB(5); “1-I coefficienti di AUTOCORRELAZIONE”
PRINT TAB(5); “2-Il test di DURBIN-WATSON, che misura la correlazione interna”
PRINT TAB(5); “3-Il nuovo test di LIN-MUSHOLKAR per la normalit…”
PRINT TAB(5); “4-L’ANALISI SPETTRALE LINEARE per il periodogramma”
PRINT : PRINT TAB(5); ” Nei punti 2 e 3 consultare le tabelle; dei’cutoff values'”
2 IF INKEY$ = “” THEN 2
CLS
DIM x(100), m(100), sd(100), L(100), rh(100), e(100), g(100), f(100)
INPUT “Immetti il numero dei dati “, n
INPUT “Immetti il numero dei lags h “, n1
‘FOR i = 1 TO n
‘PRINT “x(“; i; : INPUT “)=”, x(i): f(i) = x(i)
‘s(0) = s(0) + x(i)
‘NEXT i
‘GOTO 14

PRINT “CALCOLO DEI COEFFICIENTI DI AUTOCORRELAZIONE”
PRINT
CLS : PRINT “TABELLA DEI COEFFICIENTI DI AUTOCORRELAZIONE”: PRINT
LPRINT “TABELLA DEI COEFFICIENTI DI AUTOCORRELAZIONE”: LPRINT
FOR i = 1 TO n
READ x(i): f(i) = x(i)
s(0) = s(0) + x(i)
NEXT i
xm = s(0) / n
PRINT
PRINT “LAG h”, “COEFF. AUTOCORRELAZIONE”
PRINT
FOR h = 1 TO n1
FOR t = 1 TO n – h
so = (x(t) – xm) * (x(t + h) – xm): s(1) = s(1) + so
NEXT t
FOR t = 1 TO n
so = (x(t) – xm) ^ 2: s(2) = s(2) + so
NEXT t
rh(h) = n * s(1) / ((n – h) * s(2))
PRINT h, rh(h)
s(1) = 0: s(2) = 0
NEXT h
er = 2 / (SQR(n))
PRINT
PRINT “ERROR +/- “, er
s(0) = 0: s(1) = 0: s(2) = 0
PRINT
73 a$ = INKEY$: IF a$ = “” THEN 73
23 LPRINT “LAG H”, “COEFF.AUTOCORRELAZIONE”: LPRINT
FOR i = 1 TO n1: LPRINT USING “###”; i;
LPRINT USING ” #.###”; rh(i): NEXT i
LPRINT : LPRINT “er=”, er

GOTO 20

13 RESTORE
CLS : s(1) = 0: s(2) = 0
PRINT : PRINT “LIN-MUDHOLKAR TEST PER LA GAUSSIANA (TEST DI NORMALITA’)”
21 FOR a = 1 TO n
FOR i = 1 TO n
b = i
IF i = a THEN READ x(i): k = 1: GOTO 7
READ x(b): v = x(b)
IF k = 1 THEN b = i – 1
x(b) = v
7 NEXT i
k = 0
‘FOR i = 1 TO n – 1
‘PRINT x(i); : NEXT i
GOSUB 1000
s(1) = 0: s(2) = 0
RESTORE
NEXT a

‘SERIE STORICA ORIGINALE eliminati 3 autliers
‘ DATA .033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043
‘ DATA .066,.053,.082,.06,.08,.076,.056,.036,.05,.053,.056,.058,.061,.063,.065
‘ DATA .068,.0815,.095,.079,.063,.069,.074,.08,.0765,.073,.0695,.066,.093,.083
‘ DATA .073,.063,.074,.067,.06,.086,.08,.073,.067,.089,.064,.087,.079,.07,.065
‘ DATA .06,.063

‘EFFETTO STAGIONALE (modello additivo)
‘DATA .0022,-.0030,.0002,-.0053,.0070,.0026,.0107,.0054,-.0042,-.0083,-.0037
‘DATA -.0075,.0022,-.0030,.0002,-.0053,.0070,.0026,.0107,.0054,-.0042,-.0083
‘DATA -.0037,-.0075,.0022,-.0030,.0002,-.0053,.0070,.0026,.0107,.0054,-.0042
‘DATA -.0083,-.0037,-.0075,.0022,-.0030,.0002,-.0053,.0070,.0026,.0107,.0054
‘DATA -.0042,-.0083,-.0037,-.0075,.0022,-.0030,.0002,-.0053,.0070,.0026,.0107
‘DATA .0054,-.0042,-.0083,-.0037,-.0075
DATA 30,29,34,61,54,63,62,41,43,30,35,19,21,25,36

s(1) = 0: s(2) = 0: s(3) = 0: s(4) = 0: s(5) = 0

FOR i = 1 TO n

s(1) = s(1) + m(i): ‘Sum(X)

s(2) = s(2) + m(i) ^ 2: ‘SUM((X)^2)

s(3) = s(3) + m(i) * L(i): ‘SUM(X*Y)

s(4) = s(4) + L(i): ‘SUM(Y)

s(5) = s(5) + L(i) ^ 2: ‘SUM((Y)^2)

NEXT i
PRINT
‘PRINT s(1); s(2); s(3); s(4); s(5)
r = (n * s(3) – s(1) * s(4)) / (SQR(n * s(2) – s(1) ^ 2))
r = r * SQR(n * s(5) – s(4) ^ 2)

PRINT
PRINT “r=”, r
PRINT : PRINT
74 IF INKEY$ = “” THEN 74
LPRINT “STATISTICA DI LIN-MUDHOLKAR PER LA GAUSSIANA”
LPRINT : LPRINT “r=”, r
8 IF INKEY$ = “” THEN 8 ELSE 14
20 PRINT
PRINT “TEST STATISTICO DI DURBIN-WATSON PER L’AUTOCORRELAZIONE”
PRINT : PRINT
s(1) = 0: s(2) = 0
FOR i = 2 TO n
y = x(i) – x(i – 1)
s(1) = s(1) + y ^ 2
NEXT i
FOR i = 1 TO n
y = x(i) ^ 2
s(2) = s(2) + y
NEXT i
DW = s(1) / s(2)
PRINT “DW= “, DW
72 a$ = INKEY$: IF a$ = “” THEN 72
LPRINT : LPRINT “STATISTICA DI DURBIN-WATSON PER L’AUTOCORRELAZIONE”
LPRINT : LPRINT “DW=”, DW
LPRINT
1 a$ = INKEY$: IF a$ = “” THEN 1 ELSE 13
14 PRINT : PRINT ” ANALISI DI FOURIER E PERIODOGRAMMA”: PRINT
t = n
p1 = INT((n – 1) / 2) ‘max frequenza
c1 = COS(2 * 3.1415926# / t)
s1 = SIN(2 * 3.1415926# / t)
s = 0: p = 0
c = 1
210 u1 = 0: u2 = 0
k = t
230 u = f(k) + 2 * c * u1 – u2
u2 = u1
u1 = u
k = k – 1
IF k 0 THEN 340
a = a / 2
PRINT “COEFFICIENTI DI FOURIER”
PRINT TAB(4); “P”; TAB(12); “ALFA”; TAB(21); “BETA”
LPRINT : LPRINT “COEFFICIENTI DI FOURIER”: LPRINT
LPRINT TAB(3); “P”; TAB(12); “ALFA”; TAB(21); “BETA”: LPRINT
340 PRINT USING “##.##^^^^”; p; a; b: LPRINT USING “###”; p;
LPRINT USING ” ##.##^^^^”; a; b
IF p = 0 THEN 480
e(p) = SQR(a * a + b * b)
t2 = ABS(a / b)
t2 = 360 / 2 / 3.1415926# * ATN(t2)
IF b > 0 THEN 450
IF a > 0 THEN 430
t2 = 180 + t2
GOTO 470
430 t2 = 180 – t2
GOTO 470
450 IF a > 0 THEN 470
t2 = 360 – t2
470 g(p) = t2
480 IF p = p1 THEN 540
q = c1 * c – s1 * s
s = c1 * s + s1 * c
c = q
p = p + 1
GOTO 210
540 PRINT : PRINT
PRINT ” ANALISI ARMONICA”: PRINT
LPRINT : LPRINT “ANALISI ARMONICA”: LPRINT
PRINT ” FR”; TAB(6); “FREQ1”; TAB(15); “PERIOD”; TAB(23); “AMPIEZZA”; TAB(34); “FASE”
LPRINT ” FR”; TAB(6); “FREQ1”; TAB(15); “PERIOD”; TAB(23); “AMPIEZZA”; TAB(34); “FASE”
LPRINT
FOR i = 1 TO p1: k1 = n / i: k2 = i / n
PRINT USING “###”; i;
PRINT USING “##.##^^^^”; k2; k1; e(i); g(i)
LPRINT USING “###”; i;
LPRINT USING “##.##^^^^”; k2; k1; e(i); g(i)
NEXT i
75 IF INKEY$ = “” THEN 75
CLS
15 END

1000 FOR i = 1 TO n – 1
s(1) = s(1) + x(i)
NEXT i
m(a) = s(1) / (n – 1)
FOR i = 1 TO n – 1
y1 = x(i) – m(a)
s(2) = y1 ^ 2 + s(2)
NEXT i
sd(a) = SQR(s(2) / (n – 2))
L(a) = sd(a) ^ (2 / 3)
‘PRINT m(A); sd(A); l(A)
RETURN

SCRIPTS IN R PER IL CALCOLO DEI COEFFICIENTI DI CORRELAZIONE DELLA SERIE STORICA REALE yt IN VIA DI ANALISI A cura del dott. Piero Pistoia

QUESTO PARAGRAFO E’ in via di correzione!

# Intanto trascriviamo nel vettore yt i 60 dati della conc. As da cui partire. Impariamo poi a calcolare con R gli altri 5 vettori dati che faranno parte dell’analisi della nostra serie
reale e quindi della nostra esercitazione. Calcoliamo come primo vettore Mt (media mobile di ordine  12 su yt.

yt=c(.033,.043,.051,.059,.061,.063,.053,.036,.046,.056,.063,.048,.053,.043,.066,.053,.082,.06,.08,.076,.056,.036,.05,.053,
.056,.058,.061,.063,.065,.068,.0815,.095,.079,.063,.069,.074,.08,
.0765,.073,.0695,.066,.093,.083,.073,.063,.074,.067,.06,.086,.08,.073,.067,.089,.064,.087,.079,.07,.065,.06,.063)

t=1

#Come primo passo grafichiamo i dati e osserviamo se ci sono regolarità all’interno (trend, oscillazioni), precisiamo le ipotesi con un correlogramma ed un periodogramma, I dati sono mensili: Ipotizziamo comunque una oscillazione di periodo 12.

# Calcoliamo, come primo vettore, Mt (media mobile centrata e pesata di ordine 12 su yt).

yt=as.vector(yt) ; n=length(yt); Mt=c()
for(t in 7:n){Mt[t] = ((yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+
yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/13}

Mt # non gira! dice che esiste una } in più. Ris=0! In effetti, non so perchè, è apparsa una tonda in più, che continua ad apparire anche se corretta!

for(t in 7:n){Mt[t] = (yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+
yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+(yt[t+6])/2)/13}

mt=filter(yt,filter=rep(1/13,13))
# calcolo della Mm col comando filter di R: confrontare i due risultati
mt #OK

# in Mt ci sono i 48 (60-12) dati Media mobile di yt, da cui costruisco i 12 Fattori Stagionali (FStag) facendo la media dei 4 gennaio, dei 4 febbraio ecc. a partire da luglio, perchè Mt iniziava con luglio.
FSTag0=matrix(Mt, ncol=12, byrow=T)
# matrice di 4 righe (valori dei 12 mesi dei 4 anni) e 12 colonne con in ognuna le 4 conc. dei mesi dello stesso nome a partire da un luglio.
FStag1=colMeans(FSTag0)
#  in FStag1 trovo le 12 medie dei 4 mesi dello stesso nome (inizio luglio, fine giugno)
FStag=c( ,FStag1[7:12], FStag1[1:6]) # da controllare! Ordino da gennaio. NON OK!

Continua ad apparire una , in più!!

FStag=c( FStag1[7:12], FStag1[1:6]) # Ordino da gennaio.  OK!

FStag=c( FStag1[7:12], FStag1[1:6]) # da controllare! Ordino da gennaio. OK, ora sì
ESAs=rep(FSTAG,5) # EFFETTO STAGIONALE As #Nonostante le correzioni continua a scrivere la variabile FSTAG sbagliata!

ESAs=rep(FStag,5) #OK

ESAs
Yt1=yt-ESAs # Ciclo+Trend+Random
Yt1
Yt1s=as-vector(Yt1s) # smusso Yt1 con Mm 3*3 #Nonostante le varie correzioni appare una al posto di un .

Yt1s=as.vector(Yt1s) # OK

Yt1s=c()
for(i in 1:60){Yt1s[i]=(Yt1[i-2]+2*Yt1[i-1]+3*Yt1[i]+2*Yt1[i+1]+
Yt1[i+2])/9}
Yt1s # yt1 senza random; cioè Ciclo+Trend
RD=Yt1-Yt1s # forse si tratta solo di random: il Ciclo?

#Riportiamo in una tabella i 5 vettori dell’analisi su yt

#data <- data.frame(t,yt,ESAs,Yt1,RD)

# Facciamo i 5 correlogrammi dei vettori trovati: yt, ESAs, Yt1, Yt1s, RD
coryt=acf(yt)
coryt
corESAs=acf(ESAs)
corESAs
corYt1=acf(Yt1)
corYt1s=acf(Yt1s)
corYt1s
corRD=acf(RD)
# Interessante abbinare il correlogramma con il periodogramma.

ESEMPIO GUIDA IN R SULL’ANALISI DELL’ARSENICO: STAGIONALITA’ TRIMESTRALE

-Ho 60 dati mensili iniziali della concentrazione As delle sorgenti della Carlina (Siena). Col comando (comando in s.l.)) “matrix” organizzo per righe i 60 dati in 20 righe 3 colonne. Col comando “rowMeans” calcolo le 20 medie di riga che sono le medie trimestrali che ho chiamato ‘medietrim’ (FIG.1) che consideriamo come vettore dati da analizzare. Col comando “lm” trovo la retta di regressione sulle 20 medie trimestrali (FIG.1). La retta, pur significativa, non è molto adeguata ai dati; spiega solo il 48% della variazione. Trovo i valori predetti dalla retta con la funzione “predict”. Tolgo i 20 valori predetti dai 20 iniziali ottengo una nuova serie che è quella iniziale senza il trend rettilineo  nominata ‘detrend_trim’ (FIG.2); ancora con il comando matrix riorganizzo per riga ‘detrend_trim’ in 5 righe, 4 colonne. Nelle 4 colonne ci vanno i 4 valori trimestrali nell'”anno medio”. Infatti nella prima colonna ho i 5 valori del 1°trimestre dei 5 anni; nella seconda colonna i 5 valori del 2° trimestre e così via. Con il comando “colMeans” faccio le medie di queste 4 colonne ottenendo il Fattore stagionale trimestrale, costituito da 4 termini che estendo ai 5 anni (FIG.3 e Fig.4) con il comando ‘rep’. I 20 dati risultanti costituiscono l’Effetto stagionale. Togliamo ancora dai 20 valori originali (medietrim) l’Effetto Stagionale così ottenuto: ne risulta una nuova serie che è la serie iniziale priva dell’oscillazione stagionale.  ma con trend e residui (medieAdj_trim). Su questa serie con “lm” faccio una nuova regressione lineare con aumento dell’R-q fino a 58% (FIG.5), ottenendo come risultato della regressione ‘fitadjtrim. Potrei ottenere ora i residui sottraendo la retta della FIG.5 dai valori plottati. In effetti calcolo i residui col comando “resid” e plotto i residui (FIG.6). Sui residui applico il test di D.-Watson (forse converrebbe interpolare l’elemento 11); applico il comando acf su res (concludo: correlazione assente);  osservo infine i 4 grafici finali relativi ai residui ottenuti con plot(fitadj_trim). Conclusione: ritengo l’analisi  accettabile!

SCRIPTS IN R

library(graphics)
library(tseries)
library(stats)
library(UsingR)
library(lattice)
library(lmtest)

w=c(0.033,0.043,0.051,0.059,0.061,0.063,0.053,0.036,0.046,0.056,
0.063,0.048,0.053,0.043,0.066,0.053,0.082,0.06,0.08,0.076,0.056,
0.036,0.05,0.053,0.056,0.058,0.061,0.063,0.065,
0.068,0.0815,0.095,0.079,0.063,0.069,0.074,
0.08,0.0765,0.073,0.0695,0.066,0.093,0.083,
0.073,0.063,0.074,0.067,0.06,0.086,0.08,0.073,0.067,0.089,0.064,
0.087,0.079,0.07,0.065,0.06,.063)

par(ask=T)

par(mfrow=c(1,3))

trim=matrix(w,ncol=3,byrow=T)
medietrim=rowMeans(trim)

medietrim

# FIG.1
ts.plot(medietrim,type=”l”,main=”FIG.1″) #finchè non lo sostituisco posso usare abline

w1=c(1:20)
regtrim=lm(medietrim~w1)
abline(regtrim)

summary(regtrim)

val_pred_w=predict(regtrim) #calcolo i 20 valori predetti dalla prima regressione
length(val_pred_w)
detrend_trim=medietrim-val_pred_w

#FIG.2
plot(detrend_trim,type=”l”, main=”FIG.2″)
trim1=matrix(detrend_trim,ncol=4,byrow=T)
medietrim1=colMeans(trim1)
medietrim1_5anni=rep(medietrim1,5)

#FIG.3
plot(medietrim1_5anni,type=”l”,main=”FIG.3″)

par(mfrow=c(2,2))

#FIG.4
acf(medietrim1_5anni,main=”FIG.4″)

valAdjtrim=medietrim-medietrim1_5anni #trend_ random
fitadj_trim=lm(valAdjtrim~w1)

fitadj_trim

summary(fitadj_trim)

#FIG.5
plot(valAdjtrim,type=”l”,main=”FIG.5″)
abline(fitadj_trim)

scansione0001

#ANALISI RESIDUI
dwtest(fitadj_trim, alternative=”two.sided”)
#forse potremo interpolare l’elemento 11

#FIG.6
res=resid(fitadj_trim)
plot(res,type=”l”, main=”FIG.6″)

#FIG.7
acf(res, main=”FIG.7-res”)

scansione0002

par(mfrow=c(2,2))
#FIG.8-12
plot(fitadj_trim, main=”FIG.8-12″)

scansione0004

BLOCCO_NOTE_TRIMAs

ESTENSIONE DELLA PRECEDENTE GUIDA  A DIVERSE SERIE STORICHE (dati orari, giornalieri, settimanali, mensili…):  BREVI RIFLESSIONI. FATTORI ED EFFETTI STAGIONALI

Si parte da dati orari – Si abbiano dati orari, es.,  per un anno (n=24*365 dati). Faccio le osservazioni preliminari su questi dati (grafico, correlazione, periodogramma). Con i dati orari possiamo sbizzarrirci. Potrei con matrix organizzare n in una matrice per righe di 24 colonne e 365 righe. Con rowMeans potrei ottenere le 365 medie giornaliere su cui procedere all’analisi di dati mensili: grafico medie, trend, medie meno trend (o meno media): oscillazione+res, su cui applico ancora matrix per riga per ottenere, es., 30 colonne e 12 righe. Che cosa otterrei con rowMeans? Ottengo nella prima riga la media dei primi 30 giorni (media di gennaio), nella seconda la media di febbraio e cosi via fino alla riga 12 che sarà la media di dicembre; in definitiva rowMeans mi fornisce i dodici valori mensili che se si vuole possiamo continuare ad analizzare. Ma ipotizzo ci sia una oscillazione anche all’interno dei 24 dati orari, cioè in un giorno (Fattore stagionale orario). Per questo applichiamo invece ad n iniziale, organizzato in 24 colonne e 360 righe, colMeans. Nella prima colonna ci saranno i 360 valori relativi alle ore una, nella seconda, quelle relativi alle ore 2 … nella 24-esima i dati relativi alle ore 24 di ogni giorno. Ne consegue che colMeans calcola le 24  medie  di ogni ora del giorno di tutti i 365 giorni dell’anno in studio (media di 365 valori corrispondenti alle una, alle due…). Questo si chiama Fattore Stagionale relativo alle ore per l’anno in studio. Se volessi il Fattore stagionale mensile?  Con matrix dovrei organizzare i 360 dati in 30 righe e 12 colonne. Allora nella prima colonna andranno tutti i primi di ogni mese, nella seconda tutti i due di ogni mese, nella terza i tre di ogni mese e… nella trentesima tutti trenta del mese. Con colMeans troverei il Fattore Stagionale mensile. Per ottenere i relativi Effetti stagionali orari o mensili, ripeto con ‘rep’ rispettivamente i 24 o  i 30 valori lungo l’intero anno (365 volte per le ore e 12 volte per i mesi). In generale l’ES, se c’è, si toglie dai dati iniziali, per ottenere una serie nuova senza tale effetto (senza oscillazioni), ma contenente trend+residui. Da questa serie poi si toglie il trend e si studiano i residui per controllare se il nostro processo è sostenibile.  Potrei anche cercare il Fattore e l’Effetto stagionale settimanale (365/7= 52 settimane). La grandezza in studio cambia con i giorni della settimana? Per es., gli umani muoiono o fanno l’amore più di lunedi o di sabato? DA RIVEDERE

1-PERCHE’ RITENIAMO RILEVANTE OGGI UNA COMUNICAZIONE DIDATTICO-OPERATIVA SUL METODO DEI MINIMI QUADRATI APPLICATO ANCHE AD UN POLINOMIO TRIGONOMETRICO
2-BREVE DISCUSSIONE SULL’ ‘ARGOMENTO’ DELLA FUNZIONE SENO
3-USO DELL’ARCO-TANGENTE NEI PROGRAMMI
4-COME ‘FITTARE’ UNA COMBINAZIONE DI FUNZIONI LINEARI AI DATI DI UNA SERIE STORICA.

5-MODELLO DI REGRESSIONE LINEARE SEMPLICE: PRESENTAZIONE E ANTICIPAZIONI NECESSARIE’
6-PRIMA LINEA DI RICERCA NELLA REGRESSIONE LINEARE SEMPLICE
7-COME SI FA A VEDERE SE QUESTO MODELLO LINEARE E’ ACCETTABILE CON R. ANALISI DEI RESIDUI.
8-VARIE INFERENZE STATISTICHE CON R DOPO AVER ACCETTATO IL MODELLO CHE FITTA I DATI

9-COME SI COSTRUISCONO LE BANDE DI CONFIDENZA
10-SCRIPTS DEL PROGRAMMA IN R RELATIVO ALLA REGRESSIONE LINEARE SEMPLICE
11-COME SI FA A VEDERE SE QUESTO MODELLO E’ ACCETTABILE CON IL MATHEMATICA DI WOLFRAN Scripts di Piero Pistoia

1-PERCHE’ RITENIAMO RILEVANTE OGGI UNA COMUNICAZIONE DIDATTICO-OPERATIVA SUL METODO DEI MINIMI QUADRATI APPLICATO ANCHE AD UN POLINOMIO TRIGONOMETRICO

Uno dei problemi oggi della comunicazione culturale scientifica, sia nella Scuola sia nell’extrascolastico, è la settorializzazione continua di linguaggi sempre più evoluti, per descrivere oggetti sempre più complessi e intrisi di teoria (di “spirito” per dirlo nel linguaggio delle monadi di Leibniz). Poiché anche le scelte sociali a cui il cittadino deve essere chiamato a partecipare (se non vogliamo ritrovarci nella situazione dei nuovi analfabeti di G. Holton (1)), sono sempre più condizionate dalle relazioni tecniche degli esperti dei diversi settori, le strutture preposte alla comunicazione dovranno attivarsi per studiare il problema ed individuare i saperi che più di altri siano idonei alle scelte nell’attività sociale. “Libertà è partecipazione”, mutuando l’espressione da una canzone di Gaber, è un concetto che risuona spesso in questi ultimi tempi ora rievocato da alcuni versi di B. Brecht “Controlla il conto sei tu che lo devi pagare”, ora dal “sesto degli undici suggerimenti di un Dio” descritti da G. Conte (2) (“non accettare regole che tu stesso non ti sia dato, non adeguarti mai al pensiero della maggioranza ed alle sue mode, perché l’opinione dei più spesso non è garanzia”…).
La difficoltà incontrata nell’imparare Scienza e nell’assimilarla sembra non dipendere solo da una proposta di insegnamento più o meno efficace o gradevole, ma rimanda alla natura stessa della Scienza, che appare poco congeniale alla biologia delle mente umana (L. Wolpert, 1996 (4), parla di “scienza innaturale”). Sembra esistere cioè una discontinuità profonda fra senso comune e Scienza, con scarsa possibilità di trasferire facilmente cose della scienza al senso e all’intuito comuni e forti difficoltà ad “incarnare” i concetti scientifici. Il cervello dell’Homo sapiens, evolutosi da qualche milione a quarantamila anni fa, in interazione con un ambiente di sopravvivenza “sentito” al fine della caccia e raccolta, mediate attraverso i primi strumenti litici, la prima tecnologia che riusciva a costruire, è progetto di un genoma rimasto, fin da allora pressoché invariato, in eredità alla specie . Infatti a partire da circa 40000 anni fa il cervello umano iniziò a controllare l’evoluzione: invece di continuare a modificarsi con l’ambiente, modificava l’ambiente a sé (evoluzione da biologica a culturale di Dobzansky). Così oggi possediamo pressoché lo stesso cervello del nostro lontano antenato cacciatore raccoglitore e le pulsioni ereditate sono le responsabili delle limitate potenzialità del nostro senso comune non previsto per la comprensione scientifica, ma semmai per una tecnologia sganciata dalla scienza. L’attività scientifica può procedere solo “rompendo” con la conoscenza prescientifica ed il senso comune, come fanno gli animali sagaci di Rilke, che, sagaci appunto ed avventurosi, si sentono a disagio nella situazione pacata e di certezza in cui si trovano. Non si tratta di un prolungamento o raffinamento o ampliamento del senso comune come talora si legge, ma di qualcosa di nuovo: una specie di senso per la comprensione scientifica.. “Finchè la Scuola non fa capire quale spostamento di quadri concettuali è necessario per impadronirsi di alcuni piccoli elementari pilastri della Fisica e della Biologia, la Scuola non ha risposto alle domande a cui dovrebbe rispondere” (P. Rossi, 1996, (5)).
Una possibilità ipotetica di indebolire il senso comune si propone con un approfondimento dei processi di comprensione scientifica in un insegnamento, a livello orizzontale, intensivo e sostenuto dall’uso di una programmazione il più possibile congeniale ai processi della natura (linguaggi come il Mathematica di Wolfram) e, a livello verticale, in un insegnamento a spirale (proposto sia dal primo sia dal secondo Bruner (6)) in più riprese nel tempo. Infatti, come suggerisce la Teoria del Darwnismo Neurale del neurobiologo G. Edelman, 1998 (7), all’interno del cervello sotto la pressione delle argomentazioni collegate alla comprensione dei concetti scientifici, si possono costituire nuovi circuiti cerebrali (selezionando gruppi neurali con l’aprire o chiudere sinapsi), per cui a lungo andare quel pezzo di cultura si “incarna” (senza apportare modifiche al genoma?), innescando però processi intuitivi e creativi con arricchimento del senso comune rispetto all’argomento in studio. Si viene così a favorire quel processo automanico ed inconsapevole (Einfunlung=immedesimazione) che potrebbe catalizzare ipotesi creative, ovvero indovinare il mondo.
La rapida obsolescenza dei concetti scientifici acquisiti nella scuola, le possibile tendenze riduttive delle nuove riforme che sembrano indirizzare l’insegnamento, in maniera più o meno mediata o camuffata, verso un inserimento più efficace nelle aziende, e la necessità armai stringente nella vita sociale di partecipare in modo sempre più esperto ai progetti e alle decisioni, se non vogliamo diventare cittadini tagliati fuori dalle scelte di sopravvivenza, spingerebbero verso una riformulazione dei curricola scolastici così da includere anche a livelli più bassi di scolarizzazione saperi indispensabili per queste scelte onde innescare l’insegnamento a spirale per facilitarne il trasferimento a livello di senso comune (assimilazione). Parlo dei settori culturali che riguardano per es. l’uso della Statistica, perché con essa scopriremo il profondo e talora ambiguo significato dei “grandi” numeri esprimenti misure (3), ma specialmente dell’analisi dei dati storici, sui quali e solo su essi, sarà possibile costruire oculatamente progetti e previsioni o almeno indicare con una certa probabilità i rischi per i possibili percorsi futuri.
In questo contesto proporremo, insieme a cenni sulle regressioni lineari e multilineari, già oggi oggetto dei corsi di aggiornamento per docenti di materie scientifiche, una riformulazione informatico-operativa dell’analisi di Fourier per insegnanti dell’area scientifica, privilegiando, tramite il Matematica di Wolfram e il linguaggio di R, più che le solite dimostrazioni matematiche, un modo più intuitivo e concreto di affrontare questi concetti.

(1)-Gerald Holton “Scienza, educazione ed interesse pubblico”, il Mulino, 1990.
(2)-Giuseppe Conte “Manuale di poesia” Guanda Editore,1995. pag.26
(3)-Piero Pistoia “Esempi di analisi statistica applicata”, Didattica delle Scienze, n. 180 en. 183. La Scuola, Brescia.
(4)-L. Wolpert “La natura innaturale della Scienza”, 1996, Dedalo editore.
(5)-Paolo Rossi “Intervista”, Le Scienze, aprile, 1996.
(6)-J. Bruner “La cultura dell’Educazione”, 1997, Feltrinelli.
(7)-Gerald Edelman “Darwinismo neurale. La teoria della selezione dei gruppi neuronali”, Einaudi Editore, 1995.

>BREVE DISCUSSIONE SULL’ ‘ARGOMENTO’ DELLA FUNZIONE SENO

Si tratta di una breve riflessione sulla funzione seno e sui modi diversi di scrivere il suo argomento con esercitazione al computer (le notazioni usate nello scrivere le funzioni ed i loro argomenti sono quelle proposte dal programma Mathematica), per evidenziare l’influenza di questi modi sulla forma dell’onda e allenare così l’intuito sulle varie questioni, in particolare per gli insegnanti di Scienze.
Un’onda sinusoidale può avere l’espressione generale: yt=A*Sin[K*α+φ], dove alfa varia da 0 a 2π , la costante moltiplicativa K, come si vede, non ha dimensioni e rappresenta il numero delle onde complete in 360° e phi, la fase, rappresenta uno spostamento orizzontale del grafico dell’onda.

se α=0:
per φ=π/2, y=A ampiezza massima dell’onda che parte appunto dal valore max=A, diventando un’onda del coseno a fase zero;
per φ=0, A=0 e l’onda parte da zero;
per φ diverso da 0, caso generale, y=A Sin(φ).

Nel contesto di una serie storica l’eq. precedente acquista una forma leggermente diversa. Chiamiamo n il numero dei dati sperimentali misurati ad intervalli di tempo uguali (serie storica); esso è anche il numero degli intervalli di osservazione e quindi il periodo T della serie (T=n), immaginando che esista almeno un ciclo oscillativo completo in n dati (anche se può non esserci). Allora α/t= 2π/T; α=2πt/T = 2πt/n e yt=A*sin((k/n*t)*2*π+φ). Il simbolo K, talora detto impropriamente frequenza, rappresenta il numero dei cicli completi in n dati sperimentali ed è il numeratore della frequenza: f=K/n. Per vedere che cosa accade delle onde yt quando, per es., n=64 dati, A=1 e k varia da 1 a 3, e t varia da 0 a 64 con φ=0, far girare sul Mathematica di Wolfram (8) il programmino seguente (da aggiustare ai dati), mantenendone con attenzione la struttura e dove con Table si costruisce il vettore Yti dei dati ricavati dalla funzione:

n=64;A=1;
Yti=N[Table[A Sin[(K/n t) 2 Pi],{t,0,64}]], dove a K si sostituiscono prima 1 per ottenere yt (64 valori dell’espressione), poi 2 per ricavare yt1, un vettore dati ancora di numerosità 64, e infine 3 per yt2.

Per graficare questi tre vettori dati: yg=ListPlot[y,PlotJoined->True, GridLines->{Automatic,Automatic}], dove ad y si sostituisce in successione yt, yt1, yt2

Si deve cioè trascrivere sulla console del mathematica (vers. 2.2,  4.1…),  le seguenti linee di istruzioni con attenzione (rispettando maiuscole, parentesi e distanze):

n=64;A=1;
K=1 “Si ha un’onda completa in 64 dati T=n”
yt=N[Table[A Sin(K/n t) 2 Pi],{t,0,n}]];
yg=ListPlot[yt,PlotJoined->True,GridLines->{Automatic, Automatic}]

K=2 “Si hanno due onde complete in 64 dati: T=n/2”
yt1=N[Table[A Sin(K/n t) 2 Pi],{t,0,64},{t,0,n}]];
yg1=ListPlot[yt1,PlotJoined->True,GridLines->{Automatic, Automatic}]
ytgg1=Show[yg,yg1] “Le due onde sullo stesso piano cartesiano”

K=3 “Si hanno tre onde complete in 64 dati: T=n/3”
yt2=N[Table[A Sin(K/n t) 2 Pi],{t,0,n}]];
yg2=ListPlot[yt2,PlotJoined->True,GridLines->{Automatic, Automatic}]
ytgg1g2=Show[yg,yg1,yt2] “Le tre onde sullo stesso piano cartesiano (si provi a meccanizzare con un For)”

Proviamo ad ottenere gli stessi risultati col programma R (da elaborare)

Se infine n=64, A=2, k=1 e φ=45°:
yt3=N[Table[2 Sin[(k t 2 Pi)/n+Pi/4],{t,0,n}]]
yg=ListPlot[yt3,PlotJoined->True, GridLines->{Automatic,Automatic}]
Naturalmente ognuno può inventare gli esempi che vuole ed esercitarsi a piacere, una volta acquisita la sintassi di questo linguaggio.

Proviamo ad ottenere gli stessi risultati col programma R (da elaborare)

USO DELL’ARCO-TANGENTE NEI PROGRAMMI

Una precisazione specifica, a nostro avviso rilevante, va fatta sull’uso dell’arcotangente nei programmi, anche perché questi interventi sono rivolti ad insegnanti di Scienze in generale e comunque un insegnamento a “spirale” serve sempre. Per il calcolo delle fasi φ delle armoniche è necessario appunto applicare l’arcotangente al rapporto fra i coefficienti ak e bk di Fourier. L’ArcTan opera sulla tangente di un certo angolo alfa e dovrebbe riportare a video l’angolo di partenza secondo la convenzione standard per la misura degli angoli. In effetti, salvo per alcuni linguaggi con due funzioni ArcTan, una delle quali darebbe il giusto risultato, appare in generale un angolo compreso fra –Pi/2 e +Pi/2 (Pi = π nel linguaggio di Mathematica e pi in quello di R). Salvo il caso in cui alfa alla partenza cade nel 1° quadrante, sul risultato dell’ArcTan in generale dovremo operare alcune correzioni memorizzate nei programmi proposti che vale la pena ricordare. Vediamo come.

Se α  alla partenza cadeva nel 2° quadrante (da 90° a 180°), es.97° = 1.69297 rad, la tangente (-8.14435) è negativa (sen + , cos -) e l’ArcTan lo riporta nel 4° quadrante fra –Pi/2 e 0 (cioè: -83°= -1.44862 rad), per cui dovrò aggiungere 180° per avere il valore di partenza nel 2° quadrante (cioè: -83+180=97°).

Se alfa alla partenza era nel 3° quadrante (da 180°a 270°), es. 187°=3.26377 rad, la tangente (.122785) è positiva (sen – e cos -), l’ArcTan lo riporta fra 0 e Pi/2 nel 1° quadrante (7°=.122173 rad), dovrò così aggiungere ancora 180° per riportarlo al quadrante di origine, nel 3°.

Se alfa era nel 4°, es. 280°=4.88692 rad con tangente -5.671282, cioè negativa (sen -, cos +), l’ArcTan riporta un valore fra –Pi/2 e 0 (cioè: -80°= -1.39626 rad); dovrò così aggiungere a -80, 360 per avere i 280° di partenza.

Precisiamo infine i seguenti casi particolari. Se la tangente è zero (angolo di partenza 0 o 180°, seno=0 e coseno diverso da0 imponiamo che l’ArcTan sia zero. Se l’angolo di partenza è 90 ovvero 270° (seno +1 o –1 e coseno 0), imponiamo che l’ArcTan sia rispettivamente 90° o –90° (-90+360). Se infine ak=0 e bk=0, caso frequente nell’analisi di Fourier quando certe armoniche sono assenti nei dati, imponiamo che ArcTan sia 0°. vedere le istruzioni di R per portare phi al quadrante giusto.

4-COME ‘FITTARE’ UNA COMBINAZIONE DI FUNZIONI LINEARI AI DATI DI UNA SERIE STORICA.
Regressione lineare semplice, Algebra matriciale per regressioni anche Multilineari e matrice inversa con esempi di calcolo 
(da elaborare)

MODELLO DI REGRESSIONE LINEARE SEMPLICE: PRESENTAZIONE E ANTICIPAZIONI NECESSARIE

Il punto essenziale non è tanto quello di trattare teoricamente il metodo dei minimi quadrati (aspetto culturale abbastanza scontato), ma di prendere piena consapevolezza che tale metodo è applicabile a qualsiasi nube di punti, al limite anche omogeneamente distribuiti nel piano cartesiano, e che fornisce in ogni caso risultati! Diventa necessario e quindi obbligatorio in una ricerca seria valutare quanto questo modello sia valido di fatto (F. Anscombe,1973 ” Graphics Statistical Analysis”, American Statistician 27(1).

Secondo noi, esistono tre linee di ricerca nell’affrontare i problemi posti dai modelli di regressione.

1-        La prima linea di ricerca, minimi quadrati s.s., ci permette di capire quanto il modello ‘fitta’ bene i dati sperimentali, cioè si adatta bene ad essi (la specifica grandezza calcolata è R-quadro, quadrato del Coefficiente di Correlazione, fra l’altro (vedere dopo),  fra yi, variabile dipendente e xi indipendente.

2-     La seconda linea di ricerca, una volta soddisfatti dell’adeguatezza della retta ai dati sperimentali, ci dobbiamo mettere nelle condizioni di poter ‘misurare statisticamente ‘ anche l’adeguatezza  di quella retta sperimentale all’ipotetica retta tracciata nell’Universo di tutti i campioni possibili. In un’analisi sulla retta di regressione, è necessario che vengano rispettate una serie di ipotesi relative alle Yi della popolazione da cui proviene il campione, come 1) valori normali : cioè per ogni Xi  la distribuzione degli Yi è una gaussiana; per ogni valore cioè di Xi nell’Universo esisterà nella terza dimensione una gaussiana, le cui ascisse si trovano lungo la retta passante per Xi,Yi* parallela all’asse Y. 2) Uguali varianze (omoscedasticità: le distribuzioni gaussiane di Yi, di cui al punto 1, devono avere uguale varianza.. 3) Linearità: le medie di tutte queste distribuzioni di Yi con uguale varianza, per ogni Xi, dovranno cadere sulla retta di regressione teorica. 4) Indipendenza: tutte le osservazioni devono essere indipendenti nel senso che i valori di un dato non devono influenzare gli altri. Non ci devono essere cioè altre relazioni causali all’interno dei dati eccetto quella espressa dalla retta di regressione.

Per controllare tutte queste ipotesi necessarie al modello relative agli Xi,Yi della popolazione universo si opera a posteriori sui RESIDUI che si possono misurare essendo stime degli ‘errori veri’ (Yi-Yi*). Ciò che resta dopo aver ‘fittato’ un qualsiasi modello, si dice residuo, per ogni xi,  la differenza fra i valori yri sulla retta sperimentale (da essa predetti) ed i corrispondenti osservati o misurati della variabile dipendente yi; residuo è quello che il modello non spiega. Queste ipotesi elencate sulla variabile dipendente della popolazione universo si riflettono direttamente sugli εi  che a loro volta agiscono sui residui che dovranno avere in qualche modo lo stesso comportamento che può essere misurato. Quindi i livelli di significanza, gli intervalli di confidenza e gli altri tests sulla regressione  sono sensibili a certi tipi di violazione dei residui e non possono essere interpretati nell’usuale modo, se esistono serie discrepanze nelle assunzioni e quindi sui residui. 

In un’analisi della regressione gli εi si pensano, come già accenato, come 1) normali, 2)indipendenti, 3)con media zero e 4)varianza, sigmo-quadro, costante. Se il modello è appropriato, i residui osservati, che sono stime degli errori veri, dovrebbero avere simili caratteristiche.

I valori dei residui (resi) si stimano meglio se ogni residuo viene diviso per la stima della deviazione standard dei residui  con N-1 al denominatore del radicando, ottenendo così la serie dei residui standardizzati. I residui standardizzati con media zero  e deviazione standard 1,  sono valori positivi, se sopra la media, e valori negativi se sotto la media. Così dire che un residuo è per es., -4721, ci fornisce poca informazione; se invece la sua forma standardizzata è -3.1, ciò ci suggerisce subito non solo che il valore osservato è minore di quello predetto, ma anche che quel residuo è certamente maggiore in valore assoluto alla maggior parte dei residui, essendo più di tre deviazioni standard.

Vedere nel proseguio l’utilizzo di grafici opportuni ed altri tests (correlogramma, periodogramma, test di DURBIN_WATSON e gli svariati tests per la normalità) per valutare se la nostra curva sperimentale permette, tramite l’analisi dei residui di passare al punto 3 onde fare inferenze statistiche dal campione alla popolazione sui parametri teorici del modello, intervalli di confidenza e bande di confidenza. Da notare con attenzione che prima di aver fatto l’analisi dei residui i processi di calcolo di cui al punto 1 che rimandano alla popolazione (ottenuti come outputs di programmi al computer o altro) devono essere lasciati in sospeso! A meno che, come avviene di fatto in generale, decidiamo di procedere, senza porci problemi, a ‘testare’ le nostre ipotesi sul comportamento della popolazione, tenendo presente che le nostre conclusioni saranno affidabili o meno secondo ciò che ricaviamo dall’analisi dei residui. Infatti Anscombe nel 1973 dimostrò con i suoi campioni fittizi che due serie di dati diversi, sottoposti a regressione lineare anche se danno stessi risultati (dai coefficienti della retta alle loro inferenze sulla popolazione tramite la statistica t, all’ R-quadro e sua l’inferenza tramite la statistica F e tutte le altre conclusioni inferenziali) e non essere un modello adeguato se non sono rispettate le assunzioni sui residui standardizzati N(0, sigma^2) compresa l’indipendenza.

3- Terza linea di ricerca.    Comunque sia, se dopo l’analisi dei residui accettiamo il modello, siamo in grado, come vedremo, di fare  inferenze verso la popolazione circa tutti i parametri teorici del modello oppure accettare quelle già fatte

PRIMA LINEA DI RICERCA NELLA REGRESSIONE LINEARE SEMPLICE. Un possibile racconto.

Iniziamo misurando o utilizzando N coppie di dati (xi,yi) relativi a due variabili sperimentali, che nel piano cartesiano x,y ipotizziamo statisticamente distribuiti come una retta (ipotesi suggerita dal grafico cartesiano o da altro) per cui sia applicabile il seguente modello matematico:

Yi* = β0 + β1 * xi + εi ;   (xi,Yi*) sono le coordinate di n  punti sulla retta nella popolazione, mentre (xi,yi) sono le n coppie dei punti sperimentali.

dove xi è l’ascissa di ogni punto sperimentale e yr è l’ordinata corrispondente sulla retta di regressione, mentre xi,yi corrispondono ai punti sperimentali; εi, sconosciuti, rappresentano quanto Yi* differisce da Yi nella popolazione. Se conoscessi gli  εi troverei i coefficienti β0  β1 teorici.

Se non conosciamo gli εi, riscriviamo il modello sostituendo i corrispondenti valori stimati a e b ( o b0 e b1) che riassumono gli εi

yi = a + b * xi:   si tratta di scrivere n equazioni sostituendo le coppie di  valori conosciuti a xi e yi. Trovati a e b, ricaverò yr, “ordinate fittate”, (yr=a+bxi) ; allora yr*-Yr =εi che è la quantità residua dell’iesima osservazione; mentre resi=yr-yi -> residui. Notare la differenza fra yi e yr .

dove xi,yi sono ancora le coppie di dati sperimentali, yr sono  i valori sulla retta di regressione che ricaverò con i minimi quadrati e i valori beta0 (β0) e beta1 (β1) invece devono essere stimati e calcolati dalle coppie (xi,yi) conosciute e indicati con a e b1 o  b0 e b1, che rappresentano le nostre incognite. Per stimare beta0 e beta1, cioè calcolare le loro stime (a e b o b0 e b1)) dato che non conosciamo gli εi, è necessario individuare in qualche modo una particolare retta tracciata attraverso i punti (xi,yi) segnati nel piano cartesiano. Ma in quanti modi possiamo tracciare questa retta? In infiniti modi! E’ necessario quindi formulare un’ipotesi per individuarla. Ecco l’ipotesi: la somma delle ‘distanze’ elevate al quadrato, misurate lungo l’asse y, fra ogni punto sperimentale (tanti quanto imax) ed il corrispondente sulla retta sia un minimo (metodo dei minimi quadrati). Tali distanze sono appunto i residui (resi) che corrispondono alle stime degli εi. Dobbiamo cioè minimizzare l’espressione seguente: (vedremo poi come):

Σ(yr-yi)^2 = Σ((a + b * xi)-yi)^2 o anche  Σ(resi)^2 -> minimo

Le ipotesi iniziali su εi sono Σεi=0, cioè media(εi)=0, distribuzione gaussiana con varianza di εi=σ^2 costante e gli εi, uniformemente distribuiti (senza correlazioni interne). Se avessimo ottenuto una retta di regressione esatta, gli errori  εi ed i residui sarebbero la stessa cosa; così ci aspettiamo che i residui ci dicano qualcosa  sul valore di σ.

Potevamo anche scegliere altre ipotesi alternative, come per es., considerare la distanza a novanta gradi sempre al quadrato ecc.. Questo metodo consiste quindi nel trovare le stime di beta0 e beta1, cioè a e b, che forniscano come somma delle loro differenze elevate al quadrato un valore più piccolo possibile, cioè Σi (resi ^2)=un minimo. (da rivedere). Facendo i conti con algebra diretta e con programmi (vedere nel proseguio) otteniamo le due equazioni seguenti per il calcolo di queste stime b0 e b1 (o  a=b0; b=b1):

b1 = Σi ((xi-xm)*(yi-ym))/Σi (xi-xm)^2  se x =xi-xm e y=yi-ym  si può scrivere anche come:
b1=ΣΣ xy / ΣΣ x^2       bi è anche uguale a b1=(nΣxy-ΣxΣy)/(nΣx^2-(Σx)^2)    dove. x=xi-xm); y=yi-ym; Σx^2=Sxq=Σ(xi-xm)^2; (Σx)^2=qSx=(Σ(xi-xm))^2
ym = b0 + b1 * xm

Inseriti in qualche programma (vedere dopo) i vettori y e x (per es., i valori della matrice y (dimensioni:60*1) e x (dim.:60*2) nel programma in qbasic allegato (MULTIREG), si ottengono direttamente b0 e b1.

Ora b0 e b1 sono di certo le migliori stime per le corrispondenti grandezze nella popolazione, anche se difficilmente i numeri saranno gli stessi.

R-QUADRO O COEFFICIENTE DI DETERMINAZIONE

Esiste una misura che che indica la bontà di adattamento del modello ai dati sperimentali del campione (xi,yi) che è il Coefficiente di Determinazione indicato con R-quadro, che corrisponde poi al quadrato del Coefficiente di Correlazione lineare di Pearson fra i valori della variabile dipendente yi e la variabile indipendente xi,  o fra le yi dei dati e il corrispondente valore sulla retta di regressione yr la cui formula è:

R=Σ(x-xm)(y.ym)/(N-1)SxSy dove Sx ed Sy sono le deviazioni standard delle variabili x e y e n è la numerosità del campione.

R-quadro, come abbiamo già detto, potrebbe essere interpretato anche come il quadrato del coefficiente di correlazione  fra yi e yr, valore predetto di y dalla regressione  (questa definizione è applicabile direttamente al calcolo di esso nella multiregressione lineare)

Per approfondire il significato di R-quadro calcoliamo quale proporzione della variabilità totale della y può essere ‘spiegata’ dalla x (cioè da modello). La variabilità totale della variabile dipendente (y), cioè yi – ym può essere divisa in due componenti: variabilità spiegata dalla regressione (yri- ym) e non spiegata yi-yri=resi=ei.

yi – ym = (yri- ym) + yi-yri

Elevando a quadrato i due membri (sopprimendo il doppio prodotto che è zero) e applicando l’operatore sommatoria , avremo:

Σ(yi – ym)^2 = Σ(yi-yri)^2 + Σ(yri- ym)^2 

che si legge: la somma totale dei quadrati (TOTAL SUM OF SQUARE = TSS) è uguale alla somma dei quadrati residuali (RESIDUAL SUM OF SQUARE = RESS) più la somma dei quadrati di regressione (REGRESSION SUM OF SQUARE = REGSS). Dividendo le sommatorie a destra per i rispettivi gradi di libertà (n-p-1 e p, dove p è il numero delle variabili indipendenti, 1 nel nostro caso) si ottengono la MEAN SQUARE RESIDUAL (MSRES) e la MEAN SQUARE REGRESSION (MSREG), la cui somma è la TOTAL MEAN SQUARE (TMS)

Per calcolare quale proporzione della variabilità totale è spiegata dalla regressione, basta dividere la somma dei quadrati di regressione per la somma dei quadrati totale.

Variazione relativa spiegata = Σ(yri-ym)^2 / Σ(yi-ym)^2

Ci accorgiamo che questo rapporto è uguale a R-quadro (dimostrare). Un R-quadro, per es., di 0.44 significherà che la retta di regressione spiega il 44% della variabilità di y. Se tutte le osservazioni cadono sulla linea di regressione R-quadro è 1. Se non vi è nessuna relazione lineare fra x e y, R-quadro è zero.

 Dividendo  REGSS e RESS per i rispettivi gradi di libertà (p e n-p-1, dove p è il numero delle variabili indipendenti) si ottengono la MEAN SQUARE RESIDUAL (MSRES) e la MEAN SQUARE REGRESSION (MSREG), la cui somma è la TOTAL MEAN SQUARE (TMS).

Se le assunzioni sono rispettate e sotto le condizioni che R-quadro pop.=0 (assenza di relazione lineare nella popolazione), il rapporto fra MSREG/MSRES è distribuito come la F di Fischer con p e n-p-1 gl. Se tale rapporto è elevato (variazione spiegata > variazione residuale), riportato sulla distribuzione di Fisher, cade nella zona proibita, l’ipotesi che r-quadro pop.=0 deve essere respinta.

COME SI FA A VEDERE SE QUESTO MODELLO LINEARE E’ ACCETTABILE CON R. ANALISI DEI RESIDUI. Uso dei comandi di R.

Seguiamo il processo di costruzione del programma in R. Dobbiamo avere i dati cioè il vettore x=c(x1,x2…xn) e il corrispondente y=c(y1,y2…yn). Facciamo il plot dei dati con
plot(x,y). Da notare che y nell’esempio è chiamato yt.
Cerchiamo con R i risultati della regressione, resultreg, usando i comandi lm(y~x) o simple.lm(x,y) (per quest’ultima caricare il package UsingR)
resultreg=lm(y~x) # o
resultreg=simple.lm(x,t) # usando UsingR
Aggiungiamo al plot la retta di regressione per precisare l’idea sull’ipotesi iniziale (scelta di una regressione lineare semplice), con
abline(lm(y~x))
Nell’oggetto ‘resultreg’ ci sono contenute tutte le informazioni della regressione, in forma matriciale: oltre alle indicazioni dei coefficienti, nella prima colonna si trovano le loro stime, nella seconda i relativi errori standard (SEb0, per l’intercetta e SEb1, per la pendenza, nella quarta i valori della statistica t di Student calcolabile dal campione (b0/SEb0, b1/SEb1) che ci permette di affermare se sono accettabili o meno le stime e nella quinta le probabilità che indicano dove cade ogni stima nella distribuzione di Student. Questi valori possono essere richiamati con coef(summary(resultreg))[1,1] per b1; [2,2] per SE_pensenza ecc. Al di sotto di questa matrice a 2 righe e 4 colonne appaiono il Residual Standard Error (RSE), l’R-quadro e la F di Fisher, richiamabili con il comando cefficients$r-quadro ecc. La grandezza R-quadro (coefficiente di determinazione) e il residual standard error (SE o RSE) misurano l’adeguatezza ai dati (per es., se R^2=0.44, significa che la retta di regressione spiega il 44% della variabilità di yi). Solo dopo l’analisi dei residui valuteremo se tali stime sono accettabili.
summary(resultreg) #riassume quasi tutti i dettagli.
Per vedere parti di tali informazioni si usano i comandi res(resultreg); coef(resultreg); predict(resultreg)

b0=coef(summary(resultreg))[1,1]
b0
b1=coef(summary(resultreg))[2,1]
b1
SEb1=coef(summary(resultreg))[2,2]
SEb1
SEb0=coef(summary(resultreg))[1,2]
SEb0
summary(res) fornisce:  min. 1° Quantile, Mediana, Media, 3° Quantile, max, sui residui

Una volta conosciuto, interpretando le informazioni contenute in resultreg, che il modello è appropriato ai dati (fitta bene i dati), è necessario analizzare se sono rispettate le ipotesi iniziali sui residui  che rimandano alla  validità del modello. Successivamente vedremo come si comporteranno i parametri incogniti, stimati dai dati, nella popolazione da cui viene estratto il campione.

Nel proseguio ricalcoleremo tutte queste grandezze usando comandi di più basso livello.

COME TESTARE I RESIDUI

Per vedere se vengono rispettate le assunzioni di linearità, cioè se davvero una linea retta ‘fitta’ bene i dati, e l’omogeneità della varianza (OMOSCEDASTICITA’), si possono plottare i residui  (y) contro i valori predetti dalla regressione (x). Se si presentano chiari patterns nei grafici detti, tali assunzioni possono essere violate. In un grafico fra valori predetti standardizzati (asse x) e residui standardizzati, i residui infatti dovrebbero essere distribuiti casualmente in una banda diffusa intorno ad una linea orizzontale, che passa per lo zero. Altre configurazioni di fasce dello stesso spessore  più o meno piegate prevedono assenza di lineatità. Se invece è lo sparpagliamento dei residui ad aumentare o diminuire con la variabile x o con i valori predetti, probabilmente l’assunzione di omoscedasticità non è rispettata. Altro modo efficace per testare l’omoscedasticità dei residui è quella di plottare i residui in funzione di ogni unità di tempo (per es., residui con i mesi, se si presentano patterns a imbuto, a clessidra, a farfalla… si prevede eteroscedasticità. Per controllare infine l’indipendenza dei residui si può osservare il correlogramma e somministrare il test di Durbin-Watson, contenuti insieme ad altro nel programma in Qbasic CORR (scritto da uno degli autori) allegato. Per testare la normalità dei residui possiamo usare svariati tests (chi-quadro, Kolmogorov e Smirnov..o il più recente Lim-Mudholkar (1980), inserito nel programma CORR.

ALTRO ANCORA  SUI RESIDUI CON R (da precisare):

1) Per la  normalità  si possono usare anche istogrammi, boxplot, plots normali;  Con Normal qqplot: i residui sono normali se il grafico rimane vicino alla linea.

2) La loro uniforme distribuzione spaziale (random, assenza di correlazione interna o trends), con plots dei residui VS tempo o ordinate;
3) costanza della loro varianza ad ogni xi, con plots dei residui VS tempo, ordinate yi e valori fittati;
con Residual VS fitted: si osserva la diffusione intorno alla retta y=0 e il trend.

Il comando plot farà molto di questo per noi se gli forniamo i risultati della regressione (resultreg):

plot(resultreg) #plotta  quattro grafici separati o su un solo piano cartesiano (se esplicitiamo il comando par(mfrow(2,2)).

SEGUE LA DESCRIZIONE DI QUESTI 4 GRAFICI ……GRAFICI RELATIVI ALL’ANALISI DEI RESIDUI IN R

Il comando plot(resultreg) plotta 4 grafici separati oppure sullo stesso piano cartesiano, se esplicitiamo il comando par(mfrow(2,2)).

– Per la normalità dei residui si possono usare anche istogrammi, boxplot e plots normali; con normal qqplot i residui saranno considerati normali se il grafico rimane vicino alla linea tratteggiata (vedere fig. Normal Q-Q). Osserviamo nella curva Normal Q-Q se i residui sono aggruppati e si allontanano dalla riga tratteggiata. Plottiamo le statistiche di ordine del campione ( standardized residuals) vs i quantili da una distribuzione normale norm(mean=0, sd=1) con il comando plot(resultreg, which=2). Possiamo testare la normalità anche con lo Shapiro-Will test:

shapiro.test(residuals(resultreg))

Se il p-value fornito è < di 0.05 significance level, si respinge l’ipotesi nulla H°, che i residui siano normalmente distribuiti. Da notare che il modello di regressione è robusto rispetto all’ipotesi di normalità. Sono più importanti le assunzioni di indipendenza e varianza costante.

– L’ipotesi di indipendenza dei residui è uno dei più importanti. La dipendenza si mostra nell’autocorrelazione che può essere positiva o negativa.  Si testa la loro uniforme distribuzione spaziale (random, assenza di correlazione interna o trends), con plots dei residui vs tempo o fitted values. Valori positivi dei residui sono seguiti da valori positivi e valori negativi da valori negativi. Si presenta così  un aspetto ciclico nei residui (Fig. Residuals vs fitted).

Si fornisce  anche un test statistico, il test di Durbin-Watson, la cui statistica D è calcolata anche dal programma in Qbasic CORR allegato, insieme alla tabella e spiegazione. Con R si fa un test a due lati con l’ipotesi nulla che la correlazione non sia zero. Se il p-value è superiore a, per es., 0.05 di livello di confidenza,  si respinge l’ipotesi che non sia zero la correlazione (c’è correlazione).

library(lmtest)

dwtest(resultreg, alternative=”two.sided”)

– La costanza della loro varianza ad ogni xi si usano grafici residui vs tempo, ordinate yi e valori fittati. Con Residuals vs fitted: si osservano la diffusione intorno alla retta y=0 e il trend (vedere fig. Residuals vs fitted). La Fig. Scale Location riporta sulle ordinate la radice quadrata dei residui standardizzati vs fitted value e controlla anch’esso se la la varianza è costante. In generale si dimostra che  la varianza dei residui cambia con i residui stessi, per cui conviene dividere ogni residuo per il suo errore standard (Residual Standard Error=radice della varianza = S*sqrt(1-hii)), ottenendo gli standardized residuals. Se resi sono tutti i residui, i residui standardizzati saranno: Rsi=resi/S*sqrt(1-hii) con i da 1 a n, dove hii, chiamata leverage, verrà definita più avanti. Se |Rsi! > 2 questo residuo rimanda ad un valore di y outlier. Cancelliamo da yt questo i-esimo e rifacciamo la regressione con gli n-1 dati, ottenendo un valore fittato Yr(i) senza il residuo Di=y(i)-Yri. Per cui: Var(Di) = S^2(i)/1-hii) e ti = Di/(S(i)/1-hii))

e i ti rappresentano i residui studentizzati cancellati. ti ha una distribuzione t con n-3 gl; riportiamo questi valori su tale distribuzione per decidere se si tratta di un outlier.
Ci aspettiamo una banda costante orizzontale con i fitted value, senza sventagliate in fuori o in dentro.

Tutto questo può essere fatto fare ad R (J. Kerns 2010, Cap.11)</i)

Forniamo tramite comandi di R infine un test statistico (il Breusch Pagan test) ancora per la costanza della varianza, senza entrare nel merito (per questo vedere “Introduction to Probability and Statistics Using R” di G.J. Kerns (prima ed.,pag. 270, prec. e seg.):

library(lmtest))

bptest(resultreg)

Si respingerà l’ipotesi nulla se BP ha un p-value che risulta superiore al livello di confidenza fissato (es., 0.05).

Fig. Residual vs Leverage – Coinvolge  Outliers, Leverage,  Distanza di Cook (da elaborare)

8-VARIE INFERENZE STATISTICHE CON R DOPO AVER ACCETTATO IL MODELLO CHE FITTA I DATI

Stime sulle grandezze della popolazione

Procediamo a fare le nostre inferenze a partire dal campione senza preoccuparci per ora delle assunzioni iniziali.

– Si possono così fare stime puntuali, se stimo il valore di un parametro della popolazione a partire dal campione, che possono essere considerate le migliori ipotesi singole immediate per una grandezza della popolazione

– si possono inferire valori degli errori standard dei coefficienti della retta, immaginando di estrarre moltissimi campioni dalla popolazione e calcolare per ognuno il parametro oggetto di inferenza; se le assunzioni sono rispettate, la distribuzione del parametro sarà gaussiana con media corrispondente al valore di quel parametro nella popolazione. Seguono le formule per il calcolo della stima dell’errore standard della pendenza b1 e della stima della varianza della popolazione, costante per ogni x:

SEb1=σβ1=σ/(sqrt((n-1)*Sx^2)) e  σ^2=S^2=Σ(yi-(b0+b1*x))^2/(n-2) Da notare che SEb1 stima σβ1

SEb1=S/sqrt((Σ((xi-xm)^2; Sx=standard deviation di x

S stima σ

la cui radice quadrata è l’errore standard della stima o deviazione standard dei residui.

– Spesso basandoci su grandezze del campione, è possibile calcolare un range di valori (centrato sul valore campionario), che, con una fissata probabilità include il valore della grandezza corrispondente nella popolazione. un tale range è detto INTERVALLO DI CONFIDENZA e la stima, stima per intervallo. E’ possibile calcolare un intervallo di confidenza per i valori della popolazione: per es., un intervallo di confidenza al 95% ancora per la pendenza β1:

b1 ± tscore*SEb1   dove ± tscore sono i due valori critici della t di Student, per n-2 gl e, per es.,significanza 0.05  o Intervallo di confidenza al 95%, (in generale se gl>30 la t di Student tende ad una gaussiana per cui tscore è circa 1.96). Così β1 sarà compreso fra b1-1.96*SEb1 e B1+1.96*SEb1 e β0 sarà compreso fra b0-1.96*SEb0 e b0+1.96*SEb0.

Con R: confint(resulreg, level=.95)

Un intervallo di confidenza al 95% significa che noi estraiamo campioni ripetuti da una popolazione, sotto le stesse condizioni,  e computiamo  per ognuno l’intervallo di confidenza al 95% per la pendenza di quel campione, il 95% di questi intervalli includerebbero il valore sconosciuto della pendenza della popolazione.  Naturalmente, poichè i valori veri della popolazione non sono conosciuti, non sapremo mai se quel particolare intervallo lo contenga. Da notare che se nell’intervallo di confidenza  per la pendenza non si trova lo zero, significherà che dovremo respingere l’ipotesi nulla che la pendenza sia zero a livello di significanza osservato dello 0.05 o meno.

– Si possono testare ipotesi che un parametro abbia un determinato valore nella popolazione: per es., β1=0 (β1 pendenza della popolazione) o R-quadro nella pop. = 0: nessuna relazione lineare. La stima dell’errore standard della pendenza b1, per esempio,  può servire per testare la seguente ipotesi: il valore della pendenza nella popolazione è zero (β1=0)? Può esserlo infatti nella popolazione e non nel campione. Se nella popolazione non esiste relazione lineare (β1=0), si conosce la distribuzione della statistica t=pendenza/errore standard della pendenza, calcolata su tutti i campioni estratti dalla popolazione: si tratta di una distribuzione di Student non n-2  gl. Se il valore del t del nostro campione (tc=b1/SEb1), inserito sulle ascisse di questa distribuzione campionaria dà livelli di significanza (lascia a destra un’area di probabilità inferiore ad una soglia prefissata (0.05,0.01…), allora l’ipotesi nulla (β1=0 nella popolazione) è da respingere e ci sarà effettivamente relazione lineare fra le due variabili nella popolazione. Tale prova però non ci fornisce informazioni relative a quanto la retta di regressione spieghi i dati effettivamente (lo fa, come si è visto in precedenza, R-quadro).

– Come si è visto la somma totale dei quadrati (total sum of squares=TSS) è uguale alla somma dei quadrati residuali (residual sum of squares=RESSS) più la somma dei quadrati di regressione (regression sum of squares=REGSS).

Sotto le condizioni che R-quadro pop.=0 (assenza di relazione lineare nella popolazione), il rapporto fra MSREG/MSRES è distribuito come la F di Fischer con p e n-p-1 gl. Se tale rapporto elevato (variazione spiegata > variazione residuale), riportato sulla distribuzione di Fisher, cade nella zona proibita, l’ipotesi che r-quadro pop.=0 deve essere respinta..

– Si possono così fare controlli incrociati: usando la statistica b1/SE, se non si può respingere l’ipotesi nulla: b1=0: ma allora l’intervallo di confidenza dovrà contenere lo zero; e la F di Fisher non sarà significativa, cioè R-quadro della popolazione=0 ecc..

Inferenza circa σ

Abbiamo due parametri a e b (o b0 e b1) ìricavati dai dati attraverso conti o programmi, tre parametri da stimare nel modello che sono σ,  β0 e β1;  σ, è la deviazione standard dei termini dell’errore; se avessimo una linea di regressione ‘esatta’ i termini dell’errore ed i residui coinciderebbbero. Invece quello che è vero è che σ deve essere stimato da:

S^2=Σ(yr – yi)^2/(n-2)=Σei^2/(n-2)

dove ei sono i residui. S^2 è uno stimatore unbiased di σ^2, cioè la sua distribuzione campionaria ha come media σ^2 della popolazione universo. E’ la divisione per n-2 che lo rende corretto. La radice positiva di S^2 è chiamata SEE standard error della stima o deviazione standard dei residui. Gi errori standard della pendenza e della intercetta sono stati calcolati da R nella seconda colonna dei risultati.

LO STANDARD ERROR DI b0

σb0=σ*sqrt(1/n + (xm^2/(n-1)*(S_x)^2)

dove (S_x)^2 è la varianza campionaria della variabile indipendente

delta2=(xi-mean(xi))^2
Sx2=sum((xi-mean(xi))^2)/(n-1) varianza campionaria  di  x
Se sostituisco Sx2 in (S_x)^2, ottengo:
SEb0=S*sqrt(1/n+(mean(xi))^2/Sx2)

LO STANDARD ERROR DI b1

σb1=σ/sqrt((n-1)*(Sx)^2)

PRECISIAMO ALCUNI PASSAGGI

Inferenza circa β1

R include il test per β1=0 che controlla l’ipotesi nulla (nessuna pendenza)

Lo stimatore b1 di β1,cioè la pendenza della linea di regressione nella popolazione, è pure uno stimatore unbiased (corretto, imparziale obbiettivo).

Lo standard error di b1 è:

SEb1=S/(SQRT(Σ(xi-xm)^2)  dove S è la radice dello stimatore precedente

PER TESTARE IPOTESI

La distribuzione del valore normalizzato di b1, cioè la sua distribuzione campionaria  è una statistica t::

Non vi è alcuna relazione lineare fra x e y quando la pendenza della linea di regressione è 0.

Il test usato è:

t=b1/Sb1   la sua distribuzione statistica quando le assunzioni sono rispettate e l’ipotesi non è una relazione lineare e la distribuzione di Studente con n-2 gl (degree of fredom)

Il test usato per testare l’ipotesi che l’intercetta è zero è:

t=b1/s(b0) ha la stessa distribuzione campionaria.

Queste statistiche t e i loro livelli di significanza a due code sono riportati nella matrice dei risultati R nelle ultime due colonne.

Se vogliamo testare se β1=certo valore si utilizzano le distribuzioni campionarie seguenti:

t=(b1- β1)/SEb1   che ha una distribuzione campionaria t di Student con n-2 gradi di libertà

t=(b0-β0)/SEb0

dove SEb0=S*(sqrt(sum(xi^2/n*sum(xi-xm)^2

A questo punto è facile fare un test di ipotesi per la pendenza della regressione lineare. Per es., se l’ipotesi nulla è H0: β1=w contro l’ipotesi alternativa Ha: β1≠w; allora si può calcolare la statistica campionaria t=(b1-w)/SEb1 e trovare il valore di probabilità corrispondente dalla distribuzione-t.

Si voglia fare un test per vedere se la pendenza di -1 che prevediamo è corretta. Procediamo con R nel test statistico.

>e=resid(resultreg) # i residui del modello resultreg

>b1=(coef(resultreg))[[‘x’]] #la parte x dei coefficienti

>S=sqrt(sum(e^2/(n-2))

>SEb1=S/sum((x-mean(xi))^2 dove x è il vettore dei valori

>t=(b1-(-1))/SEb1  # cioè +1: valore della statistica campionaria

>pt(t,gl,lower.tail=F) #trova la coda a destra per questo valore di t con n-2 gl

Il valore ottenuto raddoppia se il problema richiede due lati. Se la probabilità è inf., es. 0.005 ad un dato valore si respinge l’ipotesi che la pendenza =-1.

Nel comando summary R fa da solo il test per l’ipotesi β1=0 nella colonna (pr(>|t!)) alla posizione 2,4.

Inferenza circa β0= w in R 

SE(bo)=S*sqrt(sum(x^2)/n*sum(x.mean(x))^2)))

t=(bo-w)/SEb0 # coef(resultreg)[[‘(intercept)’]]

pt(t,13,lower.tail=TRUE) # la coda inferiore <w

COME SI COSTRUISCONO LE BANDE DI CONFIDENZA

Le bande di confidenza sono strisce lungo la retta di regressione che hanno un buon effetto visivo.  La retta di regressione è usata per predire, per ogni x, il valore di yo il valore medio di y. Quanto è accurata questa previsione? prevedere per ogni x un y o la media di y porta a due bande diverse. Il valore medio di y è soggetto ad una variabilità minore della previsione del singolo y. Ambedue gli intervalli sono del tipo

b0+b1*xi ± t* SEprevisione

L’errore standard per la previsione del valore medio di y per un dato xi è:

SEprevym=S*sqrt(1/n + (xi-xmi)^2/(sum(xi-xim)^2)) dove S è la deviazione standard campionaria dei residui ei (RSE).

Se stiamo invece tentando di predire un singolo valore yt allora SEprevyt cambia anche se solo leggermente e in funzione della numerosità:

SEprevyt = S*sqrt(1 + 1/n + (xi-xmi)^2/(sum(xi-xmi)^2))

Basta costruire una tabella a tre colonne per ciascuna previsione (prevym, prevyt). Per es., per prevym: yt, SEprevymlow, SEprevytup. Nella prima colonna ci vanno gli n valori previsti yt, nella seconda gli n valori yt-SEprevym e nella terza i 60 valori yt+SEprevym. se stampiamo, sullo stesso piano cartesiano, xi,yi queste tre curve otteniamo la banda di confidenza per la previsione della media e così per l’altro caso.

Anche la funzione predict aiuta a plottare le bande.

La funzione simple.lm del Package UsingR di Venable plotterà ambedue le bande di confidenza,  chiedendolo tramite il comando show.ci=T:

simple.lim(xi,yt, show.ci=T, conf.level=0.90)

USO DI PREDICT PER PLOTTARE LE BANDE DI CONFIDENZA

Per ottenere i valori predetti, si può usare in R la funzione predict che va chiamata attraverso un data.frame con dati x sulla colonna.

predict(resultreg, data.frame(x=c(50,60))) fornisce l’yt per queste due ascisse.

Per ogni x ordinato (x=sort(x)) però vogliamo 3 valori corrispondenti all’intevallo di confidenza:

predict(resultreg,data.frame(c=sort(x)) + level=0.9, interval=”confidence”)  questo fa una tabella con tre colonne FIT  LWR  UPR: per plottare la banda inferiore usiamo la seconda colonna a cui si accede con [0,2] e con points.

plot(xi,yt)

abline(resultreg)

ci.lwr=predict(resulreg, data.frame(x=sort(x)), level=0.9, interval=”confidence”)[2]

Si aggiunge la banda con points

points(sort(xi), ci.lwr, typt=l”) # o si usa line()

Aternativamente, possiamo plottare l’altra con la funzione curve come segue:

curve(predict(resultreg, data.frame(x=x), interval=”confidence”)[,3], add=T)

Attenzione però perchè la funzione curve vuole una funzione di x non i dati! E’ difficile da interrompere.

Riassumendo proviamo queste linee in R:

Summary(resultreg)
predict(resultreg, data,frame(xi=sort(xi)), level=0.9, interval=”Confidence”)
c1.lwr= predict(resultreg, data.frame(xi=sort(xi)),level=0.9,intervall=”Confidence”)[,2]
points(sort(xi), ci.lwr, type=”l”),
ci.upr=predict(resultreg, data.frame(x=x), level=9.9, interval=”Confidence”, add=T)[,3]
points(sort(xi),ci.upr, type=”l”)

SCRIPTS DEL PROGRAMMA IN R RELATIVO ALLA REGRESSIONE LINEARE SEMPLICE: da confrontare gli outputs di R e quelli del Mathematica di Wolfram (programma quest’ultimo scritto da P. Pistoia)!

library(graphics)
library(tseries)
library(UsingR)
library(lattice)
par(ask=T)

#Il seguente vettore dati è ottenuto togliendo dai 60 dati del
#vettore dati iniziale l’ESAs (effetto stagionale arsenico), ottenendo yt con #trend + random + ciclo

yt=c(0.0308,0.0460,0.0508,0.0643,0.0549,0.0604,0.0423,
0.0306,0.0502,0.0643, 0.0667,0.0555,0.0508,0.0430,
0.0658,0.0583,0.0750,0.0574,0.0693,0.0706,0.0602,
0.0443,0.0537,0.0605,0.0534,0.0610,0.0608,0.0683,
0.0580,0.0654,0.0708,0.0896,
0.0832,0.0713,0.0727,0.0815,0.0778,0.0795,0.0728,
0.0748,0.0590,0.0904,0.0723,0.0676, 0.0672,0.0823,
0.0707,0.0675,0.0838,0.0830,0.0728,
0.0723,0.0820,0.0614,0.0763,
0.0736,0.0742,0.0733,0.0637,0.0705)
par(ask=T)
options(digits=16)

ym=mean(yt)
ym

# LA REGRESSIONE con il calcolo di b0 e b1

xi=c(1:60)

n=length(xi)
n

plot(xi,yt, type=”l”)
abline(lm(yt~xi))

x=(xi-mean(xi))
y=(yt-mean(yt))

ai=x*y
Sxy=sum(ai)
Sx=sum(x)
Sy=sum(y)
xq=x^2
Sxq=sum(xq)
yq=y^2
qSx=Sx^2
b1=(n*Sxy+Sx*Sy)/(n*Sxq-qSx)

b0=sum(yt)/n-b1*sum(xi)/n
b1;b0

Yi=b0+b1*xi

resultreg=(lm(yt~xi)) # in resultreg ci vanno i risultati sulla yt originale
summary(resultreg) # trovo la maggior parte dei risultati

coef(summary(resultreg)) # trovo otto valori del risultato sotto forma di matrice
# 4 colonne: stima, SE, t-value, pr(>t)

b0=coef(summary(resultreg))[1,1]
b0
b1=coef(summary(resultreg))[2,1]
b1
SEb1=coef(summary(resultreg))[2,2]
SEb1
SEb0=coef(summary(resultreg))[1,2]
SEb0

#t0=b0/SEb0 # SEb0
#t1=b1/SEb1 # SEb1

n1=n-2
t=b1/SEb1
pt(t,n1,lower.tail=F)

#se l’ipotesi H0 è ß1=0.05 e Ha?0.05, calcoliamo a mano il p-value
ß1=0.05
t=(b1-ß1)/SEb1
pt(t,n1,lower.tail=F)

# Inferenza su b0
ß0=200

#SEb0=sqrt(Sq*sum(xi^2)/n*sum(xi-mean(xi))^2)
#standard error di b0

SEb0
t=b0/SEb0
t
pt(t,n1,lower.tail=F)

t=(b0-ß0)/SEb0 # se per ipotesi ß0=200:
t=(b0-ß0)/SEb0
pt(t,n1,lower.tail=T) # la coda più bassa (<200)

res=resid(resultreg)
summary(res)
plot(resultreg) #aspetta per confermare cambio pagina…
par(mfrow=c(2,2))
plot(resultreg)

Sq=sum(res^2)/(length(xi)-2) # stima di sigma^2

S=sqrt(Sq) #deviazione stardard dei residui
S

CONTI CHE TORNANO

# inferenze su b1
delta2=(xi-mean(xi))^2
Sx2=sum(delta2)
SEb1=S/sqrt(Sx2)
SEb1

SEb1=S/sqrt(sum((xi-mean(xi))^2))
SEb1

t=b1/SEb1 # statistica campionaria di b1
t

# Si trova il valore p rispetto alla distribuzione t, se ß1=0

n1=length(xi)-2
pt(t,n1,lower.tail=F)

#se l’ipotesi H0 è ß1=0.05 e Ha?0.05, calcoliamo a mano il p-value

ß1=0.05
t=(b1-ß1)/SEb1
pt(t,n1,lower.tail=F)

# Inferenza su b0
ß0=200

SEb0=S*sqrt(1/n+(mean(xi))^2/Sx2)
#standard error di b0

SEb0
t=b0/SEb0

t=b0/SEb0
t

t=(b0-ß0)/SEb0 # se per ipotesi ß0=200:
t=(b0-ß0)/SEb0
pt(t,n1,lower.tail=T) # la coda più bassa (<200)

#CALCOLO DI R-QUADRO E TEST  F DI FISCHER  per R-quadro_pop=0#

# Calcolo Total Sum of Square (TSS)#
TSS=sum((yt-mean(yt))^2)
TSS

# Calcolo Regression Sum of Square (REGSS)#
REGSS=sum((Yi-mean(yt))^2)
REGSS
# Calcolo Residual Sum of Square(RESSS)#
RESSS=sum(res^2)
RESSS

# La differenza TSS-(REGSS+RESSS) deve essere 0#
TSS-RESSS-REGSS

# Calcolo Mean Square Regression (MSREG)#
p=1
MSREG=REGSS/p
MSREG
# Calcolo Mean Square Residual (MSRES)#
MSRES=RESSS/(n-p-1)
MSRES
# Calcolo Total Mean Square (TMS)
TMS=TSS/(n-1)
TMS
# Calcolo l’R-quadro che è il coefficiente di correlazione
# fra xi e yt, oppure fra yr e yt, ovvero rappresenta
# la percentuale di variazione degli yt spiegata dalla
# regressione.
Rq=REGSS/TSS
Rq
# Calcolo l’R-quadro aggiustato per la popolazione
Rqagg=1-Sq/TMS
Rqagg
# Calcolo la statistica F di Fischer
F=MSREG/MSRES
F
# Calcolo la probabilità corrispondente alla # # ascissa F. Se tale probabilità è > di 0.975#
# (livello di significanza 0.025 a due code) #
# si respinge l’ipotesi nulla.#
pf(F,p,n-p-1,lower.tail=T)

# BANDE DI CONFIDENZA

par(mfrow=c(1,1))

yt=yt*100 # ricordiamoci di pensare divise per 100 le ordinate del grafico.
# poiche la pendenza non è risolvibile nel grafico, si fa dinuovo la regressione con yt*100

resultreg=(lm(yt~xi))
summary(resultreg)# da i risultati per i punti (xi,yt*100) confrontare con quelli (xi,yt)!
plot(xi,yt) #fa uno scatterplot dei dati e vi aggiunge la retta di regressione
abline(resultreg) #si costruisce la base per aggiungere le bande.

predict(resultreg, data.frame(xi=sort(xi)), level=0.9, interval=”confidence”)

ci.lwr= predict(resultreg, data.frame(xi=sort(xi)),level=0.9,interval=”confidence”)[,2]

points(sort(xi), ci.lwr, type=”l”)

ci.upr=predict(resultreg, data.frame(xi=xi), level=.9, interval=”confidence”, add=T)[,3]

points(sort(xi),ci.upr, type=”l”)

#Usiamo un comando di più alto livello del package UsingR di Venable

resultreg=simple.lm(xi,yt)
summary(resultreg)
simple.lm(xi,yt,show.ci=T,conf.level=0.95,pred=)

 

COME SI FA A VEDERE SE QUESTO MODELLO E’ ACCETTABILE CON IL MATHEMATICA DI WOLFRAN. Scripts di Piero Pistoia

Applichiamo il modello di regressione semplice, come esempio, su tre possibili vettori dati di serie storiche mutuati da una mia ricerca su dati reali relativi a 60 concentrazioni mensili (5 anni) di arsenico nelle acque potabili delle Sorgenti Onore della Carlina (montagna a cavallo di tre provincie, Siena, Firenze, Pisa)  che integra l’acquedotto dell’Alta val di Cecina (Pi). Nei remarks iniziali è riportata in breve come sono stati costruiti questi vettori all’interno della ricerca stessa, accennando anche ai processi  serviti per il loro calcolo. I vettori verranno inseriti uno alla volta eliminando le virgolette agli estremi; vengono allegati alcuni grafici costruiti dal programma. La ricerca verrà a suo tempo inserita nel blog, come esempio di analisi statistica su dati reali. Con qualche variazione sui valori dell’asse x è possibile inserire coppie di vettori x,y qualsiasi rendendo questo uno strumento efficace per testare ogni retta in ogni piano cartesiano.  Brevi remarks sono stati abbinati anche agli svariati tests statistici condotti. Le diverse linee di programma dovranno essere riscritte sulla console del Mathematica (vers.4.x o anche la  vers. 2.x ), con molta attenzione. Buon divertimento per quelli che lo vorranno fare.

SCRIPTS DI P. PISTOIA RELATIVO ALLA REGRESSIONE SEMPLICE CON IL LINGUAGGIO DEL MATHEMATICA DI WOLFRAM

“Si hanno dati mensili, per 5 anni (60 dati), di concentrazioni As delle
sorgenti Onore della Carlina, che servono l’Alta val di Cecina; così avremo 5
concentrazioni per i 5 gennaio, 5 concentrazioni per i 5 frebbraio e così
via. Si mediano poi per ogni mese questi 5 valori per ottenere le 12 medie
seguenti a partire da gennaio.
E’ interessante notare come abbiamo ottenuto un’oscillazione in 12 mesi fattori stagionali, che potremmo estendere ai 5 anni con il comando rep, ottenendo l’effetto stagionale (ancora 60 dati) che toglieremo dai dati originali yt, al fine di avere yt1 (ciclo+trend+random). Per vedere ‘cosa c’è dentro’ gli si
applica il periodogramma.”

“yt={0.307,0.301,0.324,0.312,0.363,0.348,0.385,0.359,0.314,0.294,0.309,0.299}”
“yt=yt/5”

“I 48 dati seguenti sono stati ottenuti da 60 dati mensile, sempre dell’As,

togliendo da essi il vettore dati ottenuto da una media mobile di ordine 12 eseguita sugli stessi. In tal modo si viene a togliere dai dati originali l’oscillazione di ordine dodici (vettore asf12), scoperta con un periodogramma applicato ad essi. Allora questi 48 dati, dopo la sottrazione di asf12, conterranno ancora trend e ciclo (datitc), se l’operazoione di media avrà
eliminato i random. Su essi potremo appliare un processo che porta all’individuazione di una retta di regressione”

“>yt={0.0483,0.0527,0.0533,0.0537,0.0543,0.0550,
0.0560,0.0588,0.0609,0.0605,0.0591,0.0588,0.0591,0.0598,0.0603,0.0605,0.0602,0.0598,
0.0602,0.0611,0.0628,0.0649,0.0668,0.0685,0.0704,0.0721,0.0734,
0.0742,0.0745,0.0756,0.0767,0.0758,0.0743,0.0740,0.0744,0.0738,0.0734,0.0738,
0.0740,0.0739,0.0747,0.0745,0.0734,0.0738,0.0744,0.0743,0.0736,0.0735}”

“Il vettore asf12 (48 dati) viene utilizzato per l’elaborazione dei 12 fattori stagionali, detti AsFS, (prendendo tutti i valori di gennaio diviso 4 (media dei 4 gennaio), tutti valori di febbraio diviso 4 (media dei 4 febbraio) fino al dodicesimo fattore per dicembre. Si costruisce poi il vettore effetto stagionale (AsES) di 60 dati allineando a partire dal gennaio iniziale, questi dodici fattori stagionali ripetuti 5 volte, coprendo 60 dati. I 60 dati seguenti si ricavano sottraendo dal vettore yt (60 dati) iniziale il vettore di 60 dati dell’effetto stagionale; in tal modo si libera da yt l’effetto stagionale (AsES).In y1t che rimane ci dovrebbero essere ciclo, trend e random, Su di esso  proveremo una regressione lineare semplice”

yt={0.0308,0.0460,0.0508,0.0643,0.0549,0.0604,0.0423,0.0306,0.0502,0.0643, 0.0667,0.0555,0.0508,0.0430,0.0658,0.0583,0.0750,0.0574,0.0693,0.0706,0.0602,
0.0443,0.0537,0.0605,0.0534,0.0610,0.0608,0.0683,0.0580,0.0654,0.0708,0.0896,
0.0832,0.0713,0.0727,0.0815,0.0778,0.0795,0.0728,0.0748,0.0590,0.0904,0.0723,0.067, 0.0672,0.0823,0.0707,0.0675,0.0838,0.0830,0.0728,0.0723,0.0820,0.0614,0.0763,
0.0736,0.0742,0.0733,0.0637,0.0705}

p1 = ListPlot[yt, PlotJoined -> True,
PlotRange -> Automatic, GridLines -> {Automatic, Automatic}]

x9 = N[Table[i, {i, 1, Length[yt]}]];
c1 = x9;
For [i = 1, i < Length[yt] + 1, i++, j = i*2; c = c1; yi = yt[[i]];
c = Insert[c, yi, j]; c1 = c];
d = Partition[c, 2]
c

f[x_] := Fit[d, {1, x}, x]

z1 = Transpose[d][[1]];
z2 = Transpose[d][[2]];
yi = z2;
x1 = Min[z1];
x2 = Max[z1];
n = Length[z1]
B0 = f[x] /. x -> 0
f[0]
f1 = f[x] /. x -> 1
f[1]
B1 = f1 – B0
m[list_List] := Apply[Plus, list]/n
xv[list_List] := m[(list – xmedia)^2]
xmedia = m[z1]
ymedia = m[z2]
xvar = xv[z1]

<< Statistics`DescriptiveStatistics`
<< Statistics`ContinuousDistributions`
<< Graphics`Legend`

xmedia = Mean[z1]
xvar = Variance[z1]
yvar = Variance[z2]
yr = f[x] /. x -> z1;
RES = yi – yr;

sd = RES^2;

“Calcolo Total Sum of Square (TSS)”
TSS = Apply[Plus, (yi – ymedia)^2]

“———————————–”

“Calcolo Regression Sum of Square (REGSS)”
REGSS = Apply[Plus, (yr – ymedia)^2]

“———————————–”

“Calcolo Residual Sum Of Square (RESSS)”
RESSS = Apply[Plus, sd]

“———————————–”

“Calcolo Mean Square Regression (MSREG)”
p = 1
MSREG = REGSS/p

“———————————–”

“Calcolo Mean Square Residual (MSRES)”
MSRES = RESSS/(n – p – 1)

“———————————–”

“Calcolo Total Mean Square (TMS)”
TMS = TSS/(n – 1)

“———————————–”

“Calcolo R-quadro che è il coefficiente di”
“correlazione fra x e y, oppure fra yi e yr,”
“ovvero rappresenta la percentuale di”
“variazione degli yi spiegata dalla”
“regressione”
Rq = REGSS/TSS

“———————————–”

“Calcolo l’R quadro aggiustato per stimare il”
“corrispondente parametro della popolazione”
ESS1 = Apply[Plus, sd];
ESS2 = ESS1/(n – 2);
TSS2 = TSS/(n – 1);
Rqagg = 1 – ESS2/TSS2

“———————————-”

“Calcolo la Statistica F (MSREG/MSRES)”
” che ha la distribuzione di Fisher”
F = MSREG/MSRES

“———————————–”

“Calcolo i livelli di significanza al 95%”
“(2 code) della F di Fisher. Se la F è >”
“(coda a destra) di F95% siamo autorizzati”
“a respingere l’ipotesi nulla (Rq=0 nella”
“popolazione”
ndist = FRatioDistribution[p, n – p – 1]
F95 = Quantile[ndist, (0.95 + 1)/2] // N

“———————————–”

“Calcolo anche la probabilità (Area)”

“corrispondente all’ascissa F. Se risulta”
” > di 0.975 cade nella zona proibita”
CDF[ndist, F] // N

“———————————–”

“Calcolo l’Errore Standard della Stima o”
“Deviazione Standard dei Residui. Si tratta della”
“stima della Deviazione Standard delle distribuzioni”
“della variabile dipendente per ogni xi”
ESS = Sqrt[ESS1/(n – 2)] // N

“———————————–”

“Calcolo t-critico all’inizio della coda corrispondente”
“ad un’area di 0.975 in una distribuzione di Student a”
“n-1 gradi i libertà e ad n-2”
ndist1 = StudentTDistribution[n – 1]
ndist2 = StudentTDistribution[n – 2]
t1 = Quantile[ndist1, {0.95 + 1}/2] // N
t2 = Quantile[ndist2, {0.95 + 1}/2] // N

“———————————–”

“Calcolo l’Errore Standard di B0 (SB0) e di B1 (SB1).”
” Ricavo poi il valore campionario di tB0 (B0/SB0)”
” e tB1 (B1/SB1), per controllare se B0 e B1 nella”
” popolazione hanno valore zero”
SB0 = ESS Sqrt[1/n + xmedia^2/((n – 1) xvar)]
SB1 = ESS/Sqrt[(n – 1) xvar]
tB0 = B0/SB0
tB1 = B1/SB1

“———————————–”

“Calcolo A: area (probabilità) lasciata a sinistra”
” dell’ascissa tB0 nella sua distribuzione, una ”
“Student a n-2 GL. P=1-A è la probabilità lasciata”
” a destra di tB0 (Livello di Significanza); se ”
“questa è < di 0.025 si respinge l’ipotesi nulla”
” (B0 e/o B1 nulli nella popolazione) ”
A = CDF[ndist2, tB0]
P = 1 – A

“———————————–”

“Faccio un calcolo analogo per B1”
A1 = CDF[ndist2, tB1]
P1 = 1 – A1

“———————————–”

“Calcolo ora gli Intervalli di Confidenza al 95% di”
” B0 e B1, nei quali con la probabilità del 95% ”
“cadranno i rispettivi valori della popolazione”
B0 + t2 SB0
B0 – t2 SB0

B1 + t2 SB1
B1 – t2 SB1

“———————————–”

“Calcolo l’Errore Standard della Stima funzione”
“di (x-xmedia) nel caso di pochi dati”
k[x_] := t1 ESS Sqrt[((n + 1)/n + (x – xmedia)^2/((n – 1) xvar))]

“———————————–”

ymin = Floor[f[Floor[x1]] – k[Floor[x1]]]
ymax = Ceiling[f[Ceiling[x2]] + k[Ceiling[x2]]]

Plot[{f[x], f[x] + k[x], f[x] – k[x]},
{x, Floor[x1], Ceiling[x2]},
PlotStyle ->

{{AbsoluteThickness[0.001]},
{AbsoluteThickness[0.001], AbsoluteDashing[{5, 5}]},
{AbsoluteThickness[0.001], AbsoluteDashing[{5, 5}]}},
PlotLabel -> FontForm[Rqagg “->Rq-Aggiustato”, {“Times”, 12}],
PlotLegend -> {“Curva  Regressione”, “Limite Interv Conf 95 %”,
“Limite Interv Conf 95 %”},
LegendPosition -> {0.8, 0},
AxesOrigin -> {Automatic, Automatic},
AxesLabel -> {“Igrometro di prova”, “Igrom. di riferimento”},
GridLines -> {Automatic, Automatic},

Evaluate[Epilog -> {PointSize[0.01], Map[Point, d]}]]

matstat0001

“ANALISI DEI RESIDUI”

“Con i residui è possibile controllare se le assunzioni”
“di linearità, normalità e varianza costante, richieste”
“per l’analisi della regressione lineare, sono o no soddisfatte.”
“Basta costruire il grafico dei Residui Standardizzati sulle”
” ordinate e la variabile indipendente standardizzata sulle”
” ascisse e/o quello dei Residui Studentizzati con ancora”
“la variabile indipendente standardizzata sulle ascisse”

VIS = (z1 – xmedia)/Sqrt[xvar];

RESSTA = RES/ESS;

k1 = Flatten[k[x] /. x -> z1];

n0 = Length[k1];

RESSTU = RES/k1;

n1 = Length[VIS];

n2 = Length[RESSTA];

n3 = Length[RESSTU];

c1 = VIS;

c = VIS;

w = n + 1

For[i = 1, i < w, i++,

j = i*2; c = VIS; yi = RESSTA[[i]];

c = Insert[c, yi, j]; VIS = c]

d1 = Partition[c, 2];

Length[d1];

ListPlot[d1]

For[i = 1, i < w, i++,

j = i*2; c = c1; yi = RESSTU[[i]];

c = Insert[c, yi, j]; c1 = c]

d2 = Partition[c, 2];

Length[d2];

ListPlot[d2]

matstat0003

ALGEBRA MATRICIALE PER REGRESSIONI ANCHE MULTILINEARI E MATRICE INVERSA; CON ESEMPI DI CALCOLO (da elaborare)

COME ‘FITTARE’ UNA COMBINAZIONE DI ONDE DEL SENO AI DATI DI UNA SERIE STORICA (da elaborare)

SCRIPTS del PERIODOGRAMMA  prima parte

#In ogni caso gli scripts dei programmi presentati in R possono essere trasferiti #direttamente sulla console di R con Copia-Incolla: il programma inizierà #nell’immediato a girare costruendo risultati e grafici i cui significati sono riassunti #nei remarks.

# INIZIO ANALISI DI FOURIER IN FORMA DI FUNCTION (a cura di P. Pistoia).

t=c(1:21)

yt=100+4*sin(2*t/21*2*pi-pi/2)+3*sin(4*t/21*2*pi+0)+
6*sin(5*t/21*2*pi-1.745)

m =(n-1)/2 # perchè n dispari

PRDGRAM<- function(yt,m) {

#alfa=-pi/2 -> 270°;  alfa=-1.175 rad (cioè -100°) -> 260°

n=length(yt)
t=c(1:n)

yt=as.vector(yt)

#m =(n/2-1) # perchè n pari

#m= numero amoniche da cercare nei dati

# libray(tseries)

# VALORI DEL PARAMETRO ak
a0=c(); k=0; a0=0;
for(t in 1:n){a0=a0+yt[t]*cos(2*pi*t*k/n)}
a0

a0=a0*2/n;a0=a0/2

a0

a=c();a[1:m]=0;
for(k in 1:m) {
for(t in 1:n){
a[k]=a[k]+yt[t]*cos(2*pi*t*k/n)}}
a=2*a/n

# vALORI DEL PARAMETRO bk

b=c();b[1:m]=0;b0=0;k=0
for(k in 1:m) {
for(t in 1:n){
b[k]=b[k]+yt[t]*sin(2*pi*t*k/n)}}

a <- as.vector(a)

for(i in 1:m){
if (abs(a[i]) < 1e-10) a[i]=0 else a[i]=a[i]}
a

for(i in 1:m){
if (abs(b[i]) < 1e-10) b[i]=0 else b[i]=b[i]}
b=2*b/n
b
# aMPIEZZE

ro <- sqrt(a^2 +b^2)
ro

for(i in 1:m){
if (abs(ro[i]) < 1e-10) ro[i]=0 else ro[i]=ro[i]}
ro
# CALCOLO DELLA FASE DI OGNI ARMONICA
# RIPORTANDO IL VALORE AL QUADRANTE GIUSTO
for(i in 1:m){
f2[i] <- abs(a[i]/b[i])
f2[i] <- atan(f2[i])*180/pi}
f2 =as.vector(f2)
f2

phi <- c()
for(i in 1:m){
# f2 <- abs(a[i]/b[i]);
# f2 <- atan(f2)*180/pi;
if(b[i]>0 & a[i]>0) phi[i] = f2[i];
if(b[i]<0 & a[i]>0) phi[i] = 180-f2[i];
if(b[i]<0 & a[i]<0) phi[i] = 180+f2[i];
if(b[i]>0 & a[i]<0) phi[i] = 360-f2[i];
if(b[i]==0 & a[i]==0) phi[i] = 0;
if((b[i]<0 & a[i]>0) | a[i]==0) phi[i]=0;
if(b[i]==0 & a[i]>0) phi[i]=90;
if(b[i]==0 & a[i]<0) phi[i]=360-90
}
# PHI FASE ARMONICHE

phi=as.vector(phi)
phi
param_a <-a
param_b <-b
ampiezza <- ro
fase <- phi
data <-data.frame(param_a,param_b,ampiezza, fase)
data

par(mfrow=c(1,4))

plot(a, xlab=”Armoniche = N° osclllazioni in 21 dati”)
plot(b, xlab=”Armoniche = N° osclllazioni in 21 dati”)

r0=as.vector(ro)
plot(ro,type=”l”,main=”PERIODOGRAMMA di dati”,
xlab=”Armoniche = N° osclllazioni in 21 dati”, ylab=”ampiezza”)

# for(i in 1:10000000) i=i

# lines(phi,type=”l”)

plot(phi,type=”l”,main=”PERIODOGRAMMA di dati”,
xlab=”Armoniche = N° osclllazioni in 21 dati”, ylab=”Fase”)

}

FINE SUBROUTINE ANALISI FOURIER

Richiamo function

PRDGRAM(yt,m)           # gira la function del periodogramma di yt ed esce il seguente grafico:

Il grafico è l'output della prima parte del programma relativo al periodogramma con R

SCRIPTS del PERIODOGRAMMA seconda parte

#ESERCITAZIONI COL PROGRAMMA SCRITTO IN R da Pier Francesco Bianchi e Piero Pistoia  vers. del 16-07-19011

#Nel precedente intervento simulavamo ad hoc una serie storica
#’tabellando’ 21 dati da tre funzioni del seno con costante additiva 100,
#con ampiezze rispettivamente 4,3,6 e ‘frequenze’ nell’ordine 2/21, 4/21,
#5/21 e infine fasi -pi/2, 0, -1.745.
#Abbiamo cioè calcolati in yt 21 dati attribuendo a t valori da 1 a 21
#nell’espressione scritta nel precedente listato. Tali dati rappresentavano
#proprio la nostra serie storica da sottoporre al Periodogramma. Tramite
#il nostro programma in R calcolammo allora i valori di ampiezze e fasi
#per le prime 10 armoniche riscoprendo nei dati le oscillazioni che c’erano.
#Per esercizio continuiamo a simulare serie storiche con espressioni
#di base analoghe modificandole, aggiungendo anche un trend lineare (x*t) e
#valori random onde controllare se il Periodogramma riesce a”sentire”,
#oltre alle oscillazioni armoniche, anche il trend e la componente casuale.
#Con l’istruzione ‘#’ elimineremo secondo la necessità le linee di
#programma non utilizzate per lo scopo prefissato.

#In ogni caso gli scripts dei programmi presentati possono essere trasferiti #direttamente sulla console di R con Copia-Incolla: il programma inizierà #nell’immediato a girare costruendo risultati e grafici i cui significati sono riassunti #nei remarks.

#Proveremo ad applicare il programma su 21 dati simulati dalle
#espressioni di una retta inclinata e da una serie random estratta da
#una distribuzione gaussiana. Intanto però sceglieremo una combinazione di seni
#interessanti più adatta a proseguire l’esercitazione.

#n=21

t=c(1:n)
#yt=0.5*t
# si tratta di una iperbole
#yt=c();yt[1:t]=0
#yt <- rnorm(t,0,1)
#yt=-4+ 0.5*t + rnorm(t,0,1)
#t=c(1:21)
#yt=100+4*sin(2*t/256*2*pi-pi/2)+3*sin(4*t/256*2*pi+0)+
#6*sin(5*t/256*2*pi-1.745) #analisi yt; tenendo come base questa espressione
# con armoniche basse, ro è sulla rampa alta dell’iperbole.
#yt=100+4*sin(2*t/n*2*pi-pi/2)+3*sin(4*t/n*2*pi+0)+
#6*sin(5*t/n*2*pi-1.745) + 0.1*t # analisi ytreg
#yt=100+2*sin(2*t/n*2*pi-pi/2)+sin(4*t/n*2*pi+0)+
#3*sin(5*t/n*2*pi-1.745) + rnorm(t,0,1)*2 # analisi ytrnorm: diminuiamo le
#ampiezze e aumentiamo i random
#yt=100+4*sin(2*t/n*2*pi-pi/2)+3*sin(4*t/n*2*pi+0)+
#6*sin(5*t/n*2*pi-1.745) + 0.5*t)+(rnorm(t,0,1)-1/2)) # analisi ytregrnorm
n=240

t=c(1:n);

yt <- 4*sin(5*t/n*2*pi)+2*sin(2*pi*30*t/n)+ 3*sin(40/n*t*2*pi)+0.1*t +
rnorm(n,0,1)*2    

# questa espressione con ‘frequenze’ alte (30,40) è
#più indicata a dimostrare che il Periodogramma ‘scopre’ anche trends
#e random oltre alle oscillazioni sinusoidali.

#Ora possiamo prevedere che cosa accade se togliamo una o due di queste tre #oscillazioni,
#basta far girare il programma nei diversi casi.
#In questo contesto al termine useremo invece, per esercizio, le
#tecniche di scomposizione di una serie storica:
#proviamo a ‘destagionalizzarla’ con due o tre medie mobili
#opportune (o magari col comando filter di R) per controllare che cosa
#rimane sottraendo da yt l’oscillazione relativa alla media mobile usata (che cosa accade ai random?). Potevamo anche’detrendizzarla’ prima
#con una regressione lineare, ovvero eliminare i random con una media
#mobile 3*3 ecc..

#m =(n-1)/2 # perchè n dispari

yt=as.vector(yt)

#m =(n/2-1) # perchè n pari

#Sarebbe possibile anche ‘meccanizzare’ il precedente processo.

m=50 # numero armoniche da cercare nei dati scelte manualmente

library(tseries)

PRDGRAM(yt,m)     #RICHIAMO DELLA FUNCTION e gira il periodogramma per i nuovi dati posti in yt::

Output della seconda parte del programma esercitazione sul periodogramma con R

Le linee successive sono in attesa di correzione-eliminazione se la function funziona!

#CENNI: ‘DESTAGIONALIZZAZIONE’ DEI DATI CON ‘FILTER’ E CON MEDIE MOBILI (da correggere e precisare: AD MAIORA!)

filt12 <- filter(yt,filter=rep(1/13,13))

filt12 <- as.vector(filt12);length(filt12)

x=seq(1:240)
fit12 <- lm(filt12~x)
summary(fit12)

b0=fit12$coefficients[1]

b1=fit12$coefficients[2]

y=b0+b1*x

ts.plot(filt12)

abline(fit12)

length(filt12);length(y)

for(i in 1:10000000) i=i

filt12=filt12[7:234]
y=y[7:234]

detrend <- filt12-y

acf(detrend)
for(i in 1:10000000) i=i
ts.plot(detrend)

for(i in 1:10000000) i=i

#par(mfrow=c(1,4)) # cancelliamo ‘#’ davanti a par per vedere i grafici insieme

ts.plot(yt); for(i in 1:10000000) i

filt12 <- filter(yt,filter=rep(1/13,13))

ts.plot(filt12);for(i in 1:10000000) i

filt72 <- filter(yt,filter=rep(1/73,73))

ts.plot(filt72);for(i in 1:10000000) i

filt96 <- filter(filt96,filter=rep(1/97,97))

ts.plot(filt96);for(i in 1:10000000) i

#E’ interessante notare con tre medie mobili centrate (non pesate!) di ordine #rispettivamente 12,72,96 in successione o come vi pare, possiamo eliminare da  yt #(con yt-filtn) i tre picchi (5, 30, 40) del periodogramma ‘applicato’  alla seguente #espressione  trigonometrica che simulava i dati:

# yt <- 6*sin(5*t/n*2*pi)  + 2*sin(2*pi*30*t/n) +
# 3*sin(40/n*t*2*pi) + 0.1*t + rnorm(n,0,1)*2,

# lasciando un’iperbole ‘pulita’  a rappresentare (da controllare!)
# la retta di regressione. Da ricordare come con le medie mobili
# si eliminano anche i valori casuali.

# Che cosa accade dei picchi se cambio n?

#n=100 #n° armoniche rilevanti rimarranno ancora a freq. 5, 30 40!

#Come esercizio ‘inventiamo’ una routine che calcoli una media mobile

#di ordine k (es.,k=12)

k=12
k=k/2+1 #se pari
k= k=(k-1)/2 #se k=pari

yt=as.vector(yt) #yt è il vettore su cui si applica la media mobile
ly=length(yt)
filt=c()
for(t in 7:ly){filt[t]=
(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+yt[t-1]+yt[t]+
yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/13}
filt #Si confrontino i valori delle medie mobili calcolate con
filt12 #la definizione di media mobile e con il comando filter:
# E’ da notare che le due serie di valori differiscono perchè la
# prima serie è ottenuta da una media mobile centrata e pesata.
# Eliminando la divisione per 2 dei termini estremi avremo lo
# stesso risultato (controllare).

#RIMANE DA TROVARE UN ALGORITMO CHE AUTOMATIZZI IL
#CALCOLO DELLE MEDIE MOBILI. PROVIAMO.

Dott. Piero Pistoia – Dott. Pier Francesco Bianchi