%% ESERCIZIO I (punti 20)
% Importare le serie dei prezzi dei titoli presenti nel file prezziMIB.xlsx
% in una timetable denominata con le prime 3 lettere (senza accenti) del
% proprio cognome. Ad esempio, nel caso di Marco Riani la timetable deve
% essere denominata RIA
%
% Tramite la funzione mdpattern analizzare l'eventuale presenza di missing
% values nel file. Discutere se la presenza di missing values dipende dalle
% date in cui i prezzi dei titoli sono stati rilevati oppure se ci sono
% serie che presentano un numero di missing values molto superiori alle
% altre.
%
% Eliminare le righe che presentano valori mancanti. Dalla timetable che
% non presenta missing values, estrarre le serie il cui nome inizia con la
% lettera B in una nuova timetable. Assegnare a questa nuova timetable le
% prime 3 lettere del proprio cognome con l'aggiunta del numero 1. Ad
% esempio, nel caso di Marco Riani questa timetable deve essere denominata
% RIA1
%
% Osservazione: è preferibile scrivere il programma in maniera molto
% generale nel senso che se nel file di input ci sono 5 serie che iniziano
% con la lettera B queste vengono estratte, se nel file di input ce sono 3
% queste 3 vengono estratt...
%
% Estrarre i prezzi dei titoli che iniziano con la lettera B
% a partire dal primo di aprile 2017.
%
% Su questa timetable calcolare la serie delle medie mobili con il
% metodo del livellamento esponenziale con la funzione movavg utilizzando
% come finestra temporale 50 e 150 termini.
% Suddividere la finestra grafica in p pannelli orizzontali (dove p è il
% numero di serie che iniziano con la lettera B)
% Rappresentare, per tutti i p titoli, tramite la
% funzione plot, le serie dei prezzi e le due serie
% delle medie mobili calcolate in precedenza.
%
% Aggiungere, SOLO ALL'ULTIMO GRAFICO, la legenda.
% Discutere il segnale che le serie storiche delle due medie mobili
% mostrano in corrispondenza di inizio marzo 2020
%
%
% Passare dalle serie originarie dei prezzi alle serie dei rendimenti
% logaritmici. Effetture la conversione dalla frequenza giornaliera alla
% frequenza trimestrale utilizzando come funzione di aggregazione la media.
% Verificare l'ipotesi di normalità delle serie storiche dei rendimenti
% trimestrali. Commentare il risultato del test di normalità per ciascuna
% serie storica


% ESERCIZIO II (punti 6)
% Generare una matrice di dati di dimensione 100x5
% dalla distribuzione chi quadrato con 2 gradi di libertà.
% Denominare questa matrice con la lettera X seguita dalle prime due
% lettere del cognome in minuscolo
% (ad esempio nel mio caso la matrice è denominata Xri)
% 1) Calcolare la somma dei quadrati degli elementi di questa matrice
% utilizzando l'operatore traccia
% 2) Calcolare la matrice di correlazione e commentare i p-values del test
% di assenza di correlazione. Mostrare sia la matrice di correlazione sia
% la matrice dei p-values in formato table.


% ESERCIZIO III (punti 4).
% Calcolare in una v.c. normale con media pari a 10 e varianza pari a 15
% 1) la probabiltà di osservare valori inferiori a 8
% 2) i quantili 0.32,0.33, ..., 0.80


%% Soluzione
close all
clear
% Importare le serie dei prezzi dei titoli presenti nel file prezziMIB.xlsx
% in una timetable denominata con le prime 3 lettere (senza accenti) del
% proprio cognome. Ad esempio, nel caso di Marco Riani la timetable deve
% essere denominata RIA
RIA=readtimetable("prezziMIB.xlsx");

% Tramite la funzione mdpattern analizzare l'eventuale presenza di missing
% values nel file. Discutere se la presenza di missing values dipende dalle
% date in cui i prezzi dei titoli sono stati rilevati oppure se ci sono
% serie che presentano un numero di missing values molto superiori alle
% altre.
out=mdpattern(RIA);
% La distribuzione dei missing values è concentrata in determinate date.
% Tutte le serie hanno sicuramente la stessa probabilità di osservare
% valore mancanti.

% Eliminare le righe che presentano valori mancanti.
RIA0=rmmissing(RIA);

% Dalla timetable che
% non presenta missing values, estrarre le serie il cui nome inizia con la
% lettera B in una nuova timetable. Assegnare a questa nuova timetable le
% prime 3 lettere del proprio cognome con l'aggiunta del numero 1. Ad
% esempio, nel caso di Marco Riani questa timetable deve essere denominata
% RIA1
varn=RIA0.Properties.VariableNames;
boo=startsWith(varn,"B");
RIA1=RIA0(:,boo);

% Estrarre i prezzi a partire dal primo di aprile 2017.
boo1=RIA1.Data>=datetime('01/04/2017','InputFormat','dd/MM/yyyy');
% oppure 
% boo1=RIA1.Properties.RowTimes>=datetime('2017-04-01');

RIA2=RIA1(boo1,:);

% Su questa timetable calcolare la serie delle medie mobili con il
% metodo del livellamento esponenziale con la funzione movavg utilizzando
% come finestra temporale 50 e 150 termini.
k=50;
mme50= movavg(RIA2,'exponential',k);
k=150;
mme150= movavg(RIA2,'exponential',k);

%% Andamento dei titoli tramite funzione plot
% Suddividere la finestra grafica in p pannelli orizzontali (dove p è il
% numero di serie che iniziano con la lettera B)
% Rappresentare, per tutti i p titoli, tramite la
% funzione plot, le serie dei prezzi e le due serie
% delle medie mobili calcolate in precedenza.
%
% Aggiungere, SOLO ALL'ULTIMO GRAFICO, la legenda.
% Discutere il segnale che le serie storiche delle due medie mobili
% mostrano in corrispondenza di inizio marzo 2020
p=size(RIA2,2);
tim= RIA2.Data;
figure
for j=1:p
    subplot(p,1,j)
    plot(tim,RIA2{:,j},'b',tim,mme50{:,j},'r--',tim,mme150{:,j},'k')
    xlabel('Time')
    title(RIA2.Properties.VariableNames{j})
    if j==p
        legend({'Prezzo del titolo' 'mm esponenziali di 50 termini' 'mm esponenziali di 150 termini'},'Location','best')
    end
end
% Vicino al mese di marzo 2020 (inizio della pandemia) tutte le serie
% delle medie mobili di 50 termini tagliano dall'alto verso il basso le
% serie delle medie mobili di 150 termini (segnale di vendita)


%% Dai prezzi giornalieri ai rendimenti trimestrali

% Passare dalla serie originaria alla serie dei rendimenti logaritmici
Rt = tick2ret(RIA2,'Method','Continuous');

% Conversione dalla frequenza giornaliera alla frequenza trimestrale
Rttrim=retime(Rt,'quarterly','mean');

% Verificare l'ipotesi di normalità delle serie storiche dei rendimenti
% trimestrali. Commentare il risultato del test di normalità per ciascuna
% serie storica
Pval=zeros(p,1);
for j=1:p
    yj=Rttrim{:,j};
    [~,pvaljb]=jbtest(yj);
    Pval(j,1)=pvaljb;
end
varnames="P val test di normalità";
rownam=Rttrim.Properties.VariableNames;
PvalT=array2table(Pval,"RowNames",rownam,"VariableNames",varnames);
disp('p values del test di Normalità')
disp(PvalT)
% Non possiamo rifiutare l'ipotesi nulla di normalità per i rendimenti
% trimestrali di BUZZIUNICEM. In tutti gli altri casi si rifiuta l'ipotesi
% di normalità

%% ESERCIZIO II (punti 6)
% Generare una matrice di dati di dimensione 100x5
% dalla distribuzione chi quadrato con 2 gradi di libertà.
% Per replicabilità dei risultati fissare il seed a 123.
% Denominare questa matrice con la lettera X seguita dalle prime due
% lettere del proprio cognome in minuscolo
% (ad esempio nel caso di Marco Riani la matrice è denominata Xri).
rng(123)
Xri=chi2rnd(2,100,5);

% 1) Calcolare la somma dei quadrati degli elementi di questa matrice
% utilizzando l'operatore traccia
% v. pagina 75 del testo
sommaQ=trace(Xri'*Xri);
sommaQchk=sum(Xri.^2,'all');
assert(abs(sommaQ-sommaQchk)<1e-12,"errore di progr nel calcolo della somma dei quadrati")

% 2) Calcolare la matrice di correlazione e commentare i p-values del test
% di assenza di correlazione. Mostrare sia la matrice di correlazione sia
% la matrice dei p-values in formato table.
figure
% v. pagina 241 del testo per la chiamata a corrplot
[R,Pval]=corrplot(array2table(Xri));
disp('Matrice di correlazione in formato table')
disp(R)
disp('Matrice dei p-values in formato table')
disp(Pval)

% Le variabili sono incorrelate di conseguenza mi attendo correlazioni
% campionarie vicino a zero e p-values elevati.

%% ESERCIZIO III (punti 4).
% Calcolare in una v.c. normale con media pari a 10 e varianza pari a 15
% 1) la probabiltà di osservare valori inferiori a 8
mu=10;
sigma=sqrt(15);
disp("Pr(X<8 in N(10,15))")
disp(normcdf(8,mu,sigma))

% 2) i quantili 0.32,0.33, ..., 0.80
% Fornire l'interpretazione del quantile 0.32.
quant=0.32:0.01:0.8;
Quantili=norminv(quant,mu,sigma);
disp(Quantili)
% Quantili(1)=8.1886l è il quantile 0.32 in una  v.c. normale con media
% pari a 10 e varianza pari a 15. La probabilità in una v.c. normale con
% media pari a 10 e varianza pari a 15 di ottenere un valore inferiore a
% 8.1886l è 0.32 In altri termini 
% normcdf(Quantili(1),10,3.872983346207417)
% =0.32