/*libname tds 'c:\Documents and Settings\hraissi\Bureau\doc_cours_insa\sasrais';*/
PROC IMPORT OUT= WORK.canadian
            DATAFILE= "C:\Documents and Settings\hraissi\Bureau\Sens\Canadian.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="Feuil1$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
/*proc print data=work.canadian;
run;*/
PROC IML;
use work.canadian;
read all var{rw U} into exch;

ma=1.24;
mi=1;
pa=0.02;
n = nrow(exch);
h = do(mi,ma,pa);
numb=(ma-mi)/pa+1;
CV=J(numb,1,0);
r=J(n,numb,0);
print, h;
print, CV;

do k=1 to numb;

a=J(n,n,0);
do i=1 to n;
do j=1 to n;
a[i,j]=((4*arcos(0))**(-0.5))*exp(-0.5*((exch[i,1]-exch[j,1])/h[k])**2);
end;
end;

b=J(n,1,0);
b=a[,+];
bmat=diag(b);
bmat=inv(bmat);
b=vecdiag(bmat);

Rb=J(n,n,0);
Rb=repeat(b,1,n);

L=J(n,n,0);
L=a#Rb;

r[,k]=L*exch[,2];

num=J(n,1,0);
num=exch[,2]-r[,k];
num=num#num;

den=J(n,1,1);
den=den-vecdiag(L);
den=den#den;
denmat=diag(den);
denmat=inv(denmat);
den=vecdiag(denmat);

leav=0;
leav=t(num)*den;
leav=leav/n;

CV[k]=leav;

end;/*fin boucle k*/

/*hsilv=0.069*0.9*(1/n)**0.2;0.069 est l'écart-type de la variable explicative,
comme l'écart interquartile est égal à Q=(0.91-0.81)/1.34=0.074, 
on prend le min des deux c'est à dire l'écart-type*/
/*print ,"fenêtre en utilisant silverman (JMulTI)",,hsilv ;*/

print ,"fenêtres en utilisant la validation croisée",,CV ;

rg=1;
val=CV[1];
do i=2 to numb;
if (val>CV[i]) then rg=i ;
if (val>CV[i]) then val=CV[i];
end;/*On sélectionne la bonne fenêtre dans le vecteur CV*/

bins=400;
grid = do(min(exch[,1]),max(exch[,1]),(max(exch[,1]) - min(exch[,1]))/bins);

atilde=J(bins+1,n,0);
do i=1 to bins+1;
do j=1 to n;
atilde[i,j]=((4*arcos(0))**(-0.5))*exp(-0.5*((grid[i]-exch[j,1])/h[rg])**2);
end;
end;

btilde=J(bins+1,1,0);
btilde=atilde[,+];
bmatilde=diag(btilde);
bmatilde=inv(bmatilde);
btilde=vecdiag(bmatilde);

Rbtilde=J(bins+1,n,0);
Rbtilde=repeat(btilde,1,n);

Ltilde=J(bins+1,n,0);
Ltilde=atilde#Rbtilde;

rtilde=J(bins+1,1,0);
rtilde=Ltilde*exch[,2];

CREATE result
            var{prev aus uk};
CLOSE result ;

CREATE result2
            var{h cv};
CLOSE result2;

CREATE result3
            var{grille estim};
CLOSE result3;

/* on a créé une table pour les résultats */

prevnz = J(n,3,0);
prevnz = r[,rg]||exch;

hcv=J(numb,2,0);
hcv=t(h)||CV;

pre = J(bins+1,2,0);
pre = t(grid)||rtilde;

print prevnz;
EDIT result ;  Append from prevnz; CLOSE result ;
EDIT result2 ;  Append from hcv; CLOSE result2 ;
EDIT result3 ;  Append from pre; CLOSE result3 ;
run;
quit;
proc sort data=result out=makemoney;
by aus;
run;
  DATA datagraph ; 
      MERGE result3 makemoney;  
 run: 
proc gplot data = datagraph;
plot estim*grille uk*aus='*' /overlay ;
symbol1 color = red interpol = join;
run;
proc gplot data = result2;
plot CV*h;
symbol1 color = red interpol = join;
run;/*Bias-variance trade off*/
quit;