clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Carica il file di dati che vogliamo analizzare
% L'opzione -ascii dice a Matlab che si tratta di
% dati in formato puro testo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load -ascii dati_reg3.txt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I dati sono stati caricati da Matlab in una matrice 
% di nome dati_reg3, cioe' con lo stesso nome del file
% Estraiamo ora la variabile y e il vettore x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1 = dati_reg3(:,1);
x2 = dati_reg3(:,2);
y = dati_reg3(:,end);

% Tracciamo lo scatter plot di x contro y
% unitamente alla retta di regressione
figure(1)
plot(x1,y,'.')
l1 = lsline;
set(l1,'Color','r','LineWidth',2)
title('Y vs. X e previsione del modello lineare')
print reg3_1.wmf -dmeta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stimiamo un modello di regressione lineare con il
% comando regress e tracciamo alcuni grafici per 
% l'analisi dei residui
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
% Prepariamo la matrice dei dati
% Il modello di regressione e': y = a0+a1*x1
X = [ones(size(x1)) x1];
[b1,b1int,res1,res1int,stats1]=regress(y,X);
% Qui troviamo il coeff. R-quadro che ci fornisce una misura 
% della bonta' di adattamento del modello ai dati
stats1(1)   
% Questo e' il grafico di figura 1
subplot(2,2,1)
[xsort,idx]=sort(x1);
plot(x1,y,'.',x1,X*b1,'.')
title(['Y e previsione vs. X - R^2=',num2str(floor(100*stats1(1))/100)])
% Grafico dei residui del modello contro il valore di x
subplot(2,2,2)
plot(xsort,res1(idx))
xlabel('X')
ylabel('Residui')
title('Residui modello vs. valori di X')
% Istogramma dei residui e curva gaussiana per
% il confronto delle distribuzioni
subplot(2,2,3)
me1=mean(res1);
sd1=std(res1);
res1std=(res1-me1)/sd1;
[fr1,xc1]=hist(res1std,40);
delta1=xc1(2)-xc1(1);
fr1 = fr1/sum(fr1.*delta1);
% Trasformo le frequenze in densita' e traccio
% l'istogramma
bar(xc1,fr1,1);
hold on;
% Traccio il grafico della guassiana
xpdf=[-5:.05:5];
plot(xpdf,normpdf(xpdf),'r');
% Test di gaussianita' dei residui con il metodo di Jarque-Bera
warning off;
[h,p]=jbtest(res1);
title(['p-value JB test = ',num2str(round(1000*p)/1000)])
% Il Q-Q Plot dei residui del modello per valutarne la gaussianita'
subplot(2,2,4)
qqplot(res1)
print reg3_2.wmf -dmeta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stimiamo un modello di regressione lineare con il
% comando regress e tracciamo alcuni grafici per 
% l'analisi dei residui
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)
% Attenzione! Aggiungiamo la variabile binaria x2 alla matrice X
% In questo modo otteniamo un modello di regressione lineare
% y = a0 + a1*x1 + a2*x2 
% Essendo x2 binaria questo modello equivale a:
% y =  a0     + a1*x1      se x2=0
% y = (a0+a2) + a1*x1      se x2=1
% Quindi si tratta di due modelli lineari con diversa intercetta
X = [ones(size(x1)) x1 x2];    
[b2,b2int,res2,res2int,stats2]=regress(y,X);
% Qui troviamo il coeff. R-quadro che ci fornisce una misura 
% della bonta' di adattamento del modello ai dati
stats2(1) 
% Tracciamo lo scatter plot di x contro y
% unitamente alla previsione del modello stimato
subplot(2,2,1)
[xsort,idx]=sort(x1);
plot(x1,y,'.',x1,X*b2,'.')
title(['Y e previsione vs. X - R^2=',num2str(floor(100*stats2(1))/100)])
% Grafico dei residui del modello contro il valore di x
subplot(2,2,2)
plot(xsort,res2(idx))
xlabel('X')
ylabel('Residui')
title('Residui modello vs. valori di X')
% Istogramma dei residui e curva gaussiana per
% il confronto delle distribuzioni
subplot(2,2,3)
% Trasformo le frequenze in densita' e traccio
% l'istogramma
me2=mean(res2);
sd2=std(res2);
res2std=(res2-me2)/sd2;
[fr2,xc2]=hist(res2std,40);
delta2=xc2(2)-xc2(1);
fr2 = fr2/sum(fr2.*delta2);
bar(xc2,fr2,1);
hold on;
% Traccio il grafico della guassiana
xpdf=[-5:.05:5];
plot(xpdf,normpdf(xpdf),'r');
% Test di gaussianita' dei residui con il metodo di Jarque-Bera
[h,p]=jbtest(res2);
title(['p-value JB test = ',num2str(round(1000*p)/1000)])
% Il Q-Q Plot dei residui del modello per valutarne la gaussianita'
subplot(2,2,4)
qqplot(res2)
print reg3_3.wmf -dmeta

% Marco Sandri
% Università di Brescia
% Dipartimento Metodi Quantitativi
% info@msandri.it