%% ESERCIZIO I
% Importare le serie dei rendimenti presenti nel file missing.xlsx.
% Tramite la funzione mdpattern analizzare l'eventuale presenza di missing
% values nel file. Discutere se la distribuzione dei missing può essere
% considerata puramente casuale.
% Selezionare le due serie che presentano un numero di valori mancanti pari
% a 11 e 33. Calcolare la matrice di correlazione tra queste due serie di
% rendimenti utilizzando solo le righe che presentano valori non mancanti.
% Commentare il pvalue del test di assenza di correlazione.


%% ESERCIZIO II 
% Importare la serie del cambio euro/renminbi (valuta cinese) dal file
% BBEX3.D.CNY.EUR.BB.AC.000.csv  
%
% Importare prima i dati in formato table nella variabile Y
% Assegnare alle due colonne i nomi
% «date» e «tasso»
% e poi convertire in time table (assegnare alla timetable il nome Yt).
% Calcolare la media dei tassi di cambio dei prezzi riferiti
% al gennaio 2019 oppure 2020
%
% Convertire la serie dalla frequenza giornaliera alla frequenza
% trimestrale.
% Calcolare sulla serie trimestrale, la serie delle medie mobili con il
% metodo del livellamento esponenziale con la funzione movavg utilizzando
% come finestra temporale 10 e 20 termini. Rappresentare, tramite la
% funzione plot, la serie trimestrale del tasso di cambio e le due serie
% delle medie mobili calcolate in precedenza. Aggiungere al grafico la
% legenda. Discutere il momento in cui si è verificato un segnale
% rialzista.
%
% Convertire la serie da frequenza giornaliera a frequenza annuale.
% Fare il grafico della serie con la funzione stackedplot
% Sulla serie annuale, calcolare i rendimenti logaritmici (rt).
% Calcolare i qqplot e i boxplot di rt (inserire i due grafici in due
% pannelli della stessa finestra grafica)
% Testare la normalità dei rendimenti logaritmici tramite il test di Jarque e
% Bera. Commentare il p-value del test.


% ESERCIZIO III
% Generare una matrice di dati di dimensione 30x3 (3 variabili X, Y e Z)
% con le seguenti caratteristiche 
% Correlazione tra X ed Y diretta e molto elevata. Correlazione inversa tra
% X e Z di discreta entità (ossia intorno a -.5 --- -0.6). 
% Un valore anomalo, eccezionalmente grande, per Z. 
%
% Tramite la funzione corrplot calcolare la matrice di correlazione e fare
% il grafico della matrice dei diagrammi di dispersione evidenziando in
% rosso le correlazioni significative (al livello di significatività
% dell'uno per mille) nei diversi diagrammi di dispersione. Regredire la
% variabile Z sulle variabili X e Y. Analizzare e commentare i residui di
% questa regressione.

%% Soluzione

Y=readtable("missing.xlsx");
out=mdpattern(Y);
% La distribuzione dei missing values è completamente non casuale.
% Le 2 variabili che presentano un numero di missing pari a 11 e 33 sono
% titolo7 e titolo3.

% boo è una matrice booleana che ha la stessa dimensione della matrice Y
% che contiene vero in corrispondenza dell'elemento i,j se l'elemento i,j è
% mancante.
boo=ismissing(Y);
% Vado a contare il numero dei missing per ogni colonna.
% nummissing è un vettore riga che in corrispondenza dell'elemento j
% contiene il numero di valori mancanti della variabile j
nummissing=sum(boo,1);

% sel è un vettore booleano che contiene true (1) quando numissing è 11
% oppure nummissing è 33 
sel=nummissing==11 | nummissing==33;
% Ysel è la table che contiene le due variabili che hanno un numero di
% valori mancanti pari a 11 e 33.
Ysel=Y(:,sel);
% Verifico che effettivamente Ysel contenga le variabili titolo3 e titolo7
head(Ysel)

% In alternativa sel poteva essere trovato come
% Vado ad escludere la prima colonna e l'ultima colonna di out
nummissingCHK=out{'totPatOrMis',2:end-1};
selout=nummissingCHK==11 | nummissingCHK==33;
Yselchk=Y(:,out.Properties.VariableNames([false selout false]));
maxdiff=max(abs(Yselchk{:,[2 1]}-Ysel{:,:}),[],'all');
assert(maxdiff<1e-15,"Le due implementazioni per trovare Ysel non coincidono")

% Calcolare la matrice di correlazione tra queste due serie di
% rendimenti utilizzando solo le righe che presentano valori non mancanti.
% Commentare il pvalue del test di assenza di correlazione.
[cor,pval]=corr(Ysel{:,:},'Rows','Complete');
disp('pvalue del test di assenza di correlazione')
disp(pval)
% Il pvalue del test di assenza di correlazione tra le due variabili è
% elevato. Accetto l'ipotesi nulla di incorrelazione tra le due variabili.
% Se la distribuzione delle due variabili è normale, l'incorrelazione
% equivale all'indipendenza.

%% Soluzione
% Importare la serie del cambio euro/renminbi (valuta cinese) dal file
% BBEX3.D.CNY.EUR.BB.AC.000.csv  
%
% Importare prima i dati in formato table nella variabile Y
% Assegnare alle due colonne i nomi
% «date» e «tasso»
% e poi convertire in time table (assegnare alla timetable il nome Yt).
Y=readtable('BBEX3.D.CNY.EUR.BB.AC.000.csv','NumHeaderLines',9);
Y=Y(:,1:2);

Y.Properties.VariableNames={'date','tasso'};
Yt=table2timetable(Y(:,2),'RowTimes',Y{:,1});

%% Calcolo la media dei tassi di cambio riferiti 
% al gennaio 2019 oppure al gennaio 2020
boo=(Yt.Time.Year==2019 | Yt.Time.Year==2020) & Yt.Time.Month==1;
Ytgen2019o2020=Yt(boo,:);
disp(['La media del tasso di cambio nel gen 2019-2020 è pari a ' num2str(mean(Ytgen2019o2020.tasso,'omitnan'))])

% Convertire la serie dalla frequenza giornaliera alla frequenza
% trimestrale.
% Calcolare sulla serie trimestrale, la serie delle medie mobili con il
% metodo del livellamento esponenziale con la funzione movavg utilizzando
% come finestra temporale 10 e 20 termini. Rappresentare, tramite la
% funzione plot, la serie trimestrale del tasso di cambio e le due serie
% delle medie mobili calcolate in precedenza. Aggiungere al grafico la
% legenda. Discutere il momento in cui si è verificato un segnale
% rialzista.

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

%% Media mobile esponenziale 
% Calcolare la serie delle medie mobili con il metodo del livellamento
% esponenziale con la funzione movavg a 10 e 20 termini.
k=10;
mme10= movavg(Ytquar,'exponential',k);
k=20;
mme20= movavg(Ytquar,'exponential',k);
tim=Ytquar.Time;

%% Rappresentazione tramite la funzione plot
% la serie trimestrale del tasso di cambio e le due serie delle medie
% mobili calcolate in precedenza.
% Inizializzo una nuova finestra grafica
figure
plot(tim,Ytquar{:,1},'b',tim,mme10{:,1},'r--',tim,mme20{:,1},'k')
legend({'Tax cambio' 'mm esponenziali di 10 termini' 'mm esponenziali di 20 termini'},'Location','best')
xlabel('Time')
%  01/01/2003 FORTE SEGNALE RIALZISTA. La media di breve (10 termini) taglia dal basso
%  all'alto la media di medio-lungo (20 termini)

%% Conversione della frequenza giornaliera alla frequenza annuale
Ytyear=retime(Yt,'yearly','mean');
figure
stackedplot(Ytyear)

%% Passare dalla serie originaria alla serie dei rendimenti logaritmici
rt = tick2ret(Ytyear,'Method','Continuous');
rt.Properties.VariableNames={'rt'};
figure
stackedplot(rt)
rtarray=rt{:,1};

%% qqplot e boxplot
figure
%  Calcolare i qqplot e i boxplot di rt (inserire i due grafici in due
% pannelli della stessa finestra grafica)
subplot(1,2,1)
% Confronto dei quantili osservati con quelli della distribuzione normale
% Se i dati sono distribuiti in maniera normale i punti si dispongono lungo
% una retta
qqplot(rtarray)
subplot(1,2,2)
boxplot(rtarray)

[h,pval]=jbtest(rtarray);
% Il p-value del test di normalità è superiore a 0.05. 
% Non posso rifiutare l'ipotesi nulla che i rendimenti annuali siano
% distribuiti in maniera normale. Occorre osservare che la lunghezza della
% serie storica è molto breve (campione molto piccolo), di conseguenza è
% meglio affermare "non posso rifiutare" piuttosto che 
% "questi rendimenti sono distribuiti in maniera normale".

% Dato il valore del MAR  (minimum acceptable return) basato
% su tasso relativo ad un titolo risk free pari a 0.004, si calcolino: 
% 1) il numero di r_t che cadono sotto la soglia del MAR (momento di ordine
% 0 oppure Shortfall risk).
% 2) La varianza dei rendimenti che cadono sotto la soglia del MAR (momento
% di ordine 2).
% 
% Definisco la soglia del MAR
MAR=0.004;

% Calcolo dei momenti di ordine 0 e 2 tramite la funzione lpm
Momp=lpm(rtarray,MAR,[0 2]);


%% ESERCIZIO III
% Generare una matrice di dati di dimensione 30x3 (3 variabili X, Y e Z)
% con le seguenti caratteristiche. 
% Correlazione tra X ed Y diretta e molto elevata. Correlazione inversa tra
% X e Z di discreta entità (ossia intorno a -.5 --- -0.6). 
% Un valore anomalo, eccezionalmente grande, per Z.
close all
rng(100)
% per replicabilità dei risultati (qualsiasi altro seed sarebbe andato
% bene)
n=30;
% Correlazione diretta molto elevata tra x e y
% Genero x in base ad una sequenza anche se potevo generarlo tramite rand
% oppure randn. 
x=(1:n)';
% Tanto più aumento coef1 tanto più l'importanza della componente casuale
% aumenta e tanto più la correlazione tra x e y diminuisce
coef1=10;
y=2+3*x+coef1*randn(n,1);
% Correlazione inversa tra X e Z di discreta entità
% Tanto più aumento coef2 tanto più la correlazione tra X e Z diminuisce
% (in valore assoluto)
coef2=100;
z=2-5*x+coef2*randn(n,1);
% Un valore anomalo (eccezionalmente grande) per Z
% Contamino un'unità a caso (in questo caso l'elemento 3 di z).
z(3)=max(z)+300;
X=[x y z];
% Tramite la funzione corrplot calcolare la matrice di correlazione
% e fare il grafico della matrice dei diagrammi di dispersione
% evidenziando in rosso le correlazioni significative (al livello di
% significatività dell'uno per mille) nei diversi diagrammi di dispersione.
Rgen=corrplot(X,'testR','on','alpha',0.001);
disp(Rgen)


% Regredisco z su x e y 
% Con il codice che segue le tre variabili della table sono denominate x, y
% e z
XX=table(x,y,z);
% L'ultima variabile della table viene presa come dipendente nella funzione
% fitlm
outLM=fitlm(XX);
disp(outLM)
% Grafico dei residui
% Considero i residui denominati studentizzati (colonna 3 della table
% outLM.Residuals) ma potevo prendere anche altre colonne
resindexplot(outLM.Residuals{:,3})
% Tutti i residui sono dentro le bande ad eccezione del valore anomalo
% (unità 3).