/* ================================================================== */
/* =========== Étape 5:  Commandes SAS pour modéliser =============== */
/* =========== avec PROC ARIMA une série chronologique ============== */
/* ================================================================== */

/* 
Remarque: Dans cette partie, vous aurez peut-être à changer que les noms de fichiers, 
          et bien entendu les modèles que vous désirez considérer.  
          Un autre changement important est le nombre associé au LEAD correspondant 
          au nombre de donnees tronquées, qui est une information utile afin d'évaluer 
          la performance prévisionnelle de votre modèle.
*/

ODS PDF FILE='ESTIMATION.pdf';
OPTIONS PS=66 LS=80 MISSING = 'N' FORMDLIM= '=';

DATA serie;
  NBOBS= _N_;
  INFILE 'serie.txt';
  INPUT Zt @@;             /* Lecture des données dans le fichier UNIX */
  Ztt = Zt;                /* Transformation des données; ici par défaut on ne transforme pas. */
;

/*============================ PROC ARIMA ===========================*/
PROC ARIMA DATA=serie;

  /*===== ETAPE 1: IDENTIFICATION =====*/
  IDENTIFY VAR=Ztt NLAG=24;             /* Étape d'identification */

 /*=====Autres exemples=====*/
 /* 
   IDENTIFY VAR=Ztt(l,12);     (1-B)(1-B**12)Zt 
   IDENTIFY VAR=Ztt(l) ;       (1-B)Zt 
   IDENTIFY VAR=Ztt(12);       (1-B**12)Zt 
   IDENTIFY VAR=Ztt(1,1);      {(1-B)**2}Zt = (1-B)(1-B)Zt = (1-2*B+B**2)Zt
*/

 /*===== ETAPE 2: ESTIMATION  =====*/
   ESTIMATE p=1 METHOD=ML PLOT PRINTALL;  /*  Étape d'estimation par
                                              la méthode du maximum 
                                              de vraisemblance */
 /*  D'autres méthodes: CLS (moindres carrés conditionnels);
                        ULS (moindres carrés inconditionnels);
 */

/*=====D'autres exemples de modeles=====*/
/*
   ESTIMATE p=3 q=(1) (4) METHOD=ML PLOT PRINTALL ;  (1- phi1*B - phi2*B**2 - phi3*B**3)Zt = (1-theta4* B**4)at;
   ESTIMATE p=(1) (3) q=2 METHOD=ULS PLOT PRINTALL;  (1-phi1*B)(1-phi3*B**3)Zt = (1-theta1B-theta2B**2)at;
   ESTIMATE p=(3) (12) METHOD=CLS PLOT PRINTALL;     (1-phi3*B**3)((1-phi12*B**12)Zt = at;
   ESTIMATE p=12 NOCONSTANT METHOD=ML ;              AR(12);
   ESTIMATE p=(12) NOCONSTANT METHOD=ML ;            (1-phi12*B**12)Zt = at.
*/

/*=====CALCUL DES PREVISIONS=====*/
   FORECAST OUT=out BACK=0 LEAD=5 ID=NBOBS; 	/* Exemple de prevision; Les résultats vont dans le DATA out */
	
/* 
  Remarque: le nombre associé à l'énoncé LEAD doit etre êgal au nombre de données dans
	      le fichier nouvelles.dat.  Ici cette valeur vaut 5.  Si vous tronquez une
            année entière de données mensuelles cela pourrait valoir 12.
*/

RUN;

/* ================================================================== */
/* ============= Etape 6: Lecture des nouvelles données  ============ */
/* ================================================================== */
/*
  Dans cette étape, vous avez formulé un modèle que vous pensez satisfaisant, 
  maintenant vous désirez en étudier la performance prévisionelle.  
  Vous allez pour ce faire effectuer des prévisions et comparer avec les 
  données qui se trouvent dans le fichier "nouvelles.txt".
*/

DATA combiner;
	MERGE serie out;
RUN;

DATA nouvelle;
  INFILE 'nouvelles.txt';
  INPUT Zt @@;
  NBOBS = _N_;
RUN;

/* =========================================================== */
/* ============= Etape 7: Premier graphique  ================= */
/* =========================================================== */
DATA retrans;
  SET combiner;
  L95orig = L95;
  U95orig = U95; 
  PREVorig = forecast;
RUN;

/*
 Remarque: Si vous avez transformé les données, les commandes du bloc précédent
           sont inexactes.
 Exemple d'utilisation: Si Zttrans = log(serie), alors:
                        L95orig = EXP(L95);
                        U95orig = EXP(U95);
                        PREVorig = EXP(forecast);
*/


DATA prevout;
  SET retrans;
  IF NBOBS <= 96;  /*  Cette valeur est le nombre d'observations de la série originale */
RUN;

GOPTIONS RESET=ALL GACCESS='sasgastd > bandes_prev.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V=DOT H=.1 C=BLACK I=JOIN L=1; 
SYMBOL2 V=DOT H=.1 C=BLACK I=JOIN L=3;
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=prevout;
  PLOT Zt*NBOBS=1 l95orig*NBOBS=2 u95orig*NBOBS=2 
       /overlay haxis=axis1 hminor=0 vaxis=axis2 vminor=0;
TITLE "SERIE ET BANDES DE PREVISIONS (ECHELLE ORIGINALE)";
run; 


/* =========================================================== */
/* ============= Etape 8: Second graphique  ================== */
/* =========================================================== */
DATA perfout;
  SET retrans;
  IF NBOBS >= 97; /*  Cette valeur devrait être n+1 où n est le nombre d'observations de la série originale */
RUN;

DATA perfout2;
  MERGE nouvelle perfout(DROP= NBOBS Zt Ztt L95 U95 FORECAST); 
RUN;

GOPTIONS RESET=ALL GACCESS='sasgastd > prev_obsfutures.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V='P' H=1 C=BLACK I=JOIN; 
SYMBOL2 V='R' H=1 C=BLACK I=JOIN;
SYMBOL3 V=DOT H=0.1 C=BLACK I=JOIN L=3;
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=perfout2;
  PLOT prevorig*NBOBS=1 Zt*NBOBS=2 l95orig*NBOBS=3 u95orig*NBOBS=3 
       /overlay haxis=axis1 hminor=0 vaxis=axis2 vminor=0;
TITLE "PREVISIONS ET BANDES DE PREVISIONS (ECHELLE ORIGINALE)";
run; 


/* ============================================================ */
/* ========= Etape 9: Performance prévisionnelle  ============= */
/* ============================================================ */
PROC PRINT DATA=perfout2;
RUN;

ODS PDF CLOSE;