/* ================================================================= */
/* ======== Méthode de Box-Jenkins : Identification ================ */
/* ================================================================= */

ODS PDF FILE='IDENTIFICATION.pdf';
OPTIONS PS=66 LS=80 FORMDLIM= '=';

/*  
Supposons que vous disposez de deux fichiers de données:
1. "serie.txt" est un fichier qui contient les données que vous désirez modéliser.
2. "nouvelles.txt" est un fichier qui contient les nouvelles données.
Ainsi, ayant vos données téléchargées d'une base de données, votre 
première étape consiste à créer ces deux fichiers, et les noms des fichiers sont
importants (si vous optez pour des noms différents vous devrez changer le code aux 
endroits pertinents).

Remarque : Le nombre de données tronquées est une information importante.  
Cette quantité doit être inscrite dans PROC ARIMA de SAS lors du calcul des prévisions.  
Voir l'énoncé LEAD.
*/

/* ================================================================== */
/* ================== Étape 1: Lecture des données ================== */
/* ================================================================== */

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

GOPTIONS RESET=ALL GACCESS='sasgastd > fig_serie.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V=DOT H=.1 C=BLACK I=JOIN L=1; 
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=serie;
  PLOT Zt*NBOBS=1 / HAXIS=axis1 HMINOR=0 VAXIS=axis2 VMINOR=0;
TITLE "SERIE ORIGINALE";
RUN; 

/* ================================================================== */
/* ========== Étape 2: Besoin de transformer les données? =========== */
/* ================================================================== */

/*
On peut utiliser comme analyse préliminaire les transformations de type racine,
inverse ou logarithmique afin de constater si la série semble posséder une variance
plus ou moins constante.
*/

DATA transfo;
  SET serie;
  ZtINV = 1/Zt;
  ZtSQRT = SQRT(Zt);
  logZt = LOG(Zt);
RUN;

GOPTIONS RESET=ALL GACCESS='sasgastd > fig_inv.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V=DOT H=.1 C=BLACK I=JOIN L=1; 
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=transfo;
  PLOT ZtINV*NBOBS=1 / haxis=axis1 hminor=0 vaxis=axis2 vminor=0;
TITLE "TRANSFORMATION INVERSE";
run; 

GOPTIONS RESET=ALL GACCESS='sasgastd > fig_sqrt.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V=DOT H=.1 C=BLACK I=JOIN L=1; 
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=transfo;
  PLOT ZtSQRT*NBOBS=1 / haxis=axis1 hminor=0 vaxis=axis2 vminor=0;
TITLE "TRANSFORMATION RACINE";
run; 

GOPTIONS RESET=ALL GACCESS='sasgastd > fig_log.pdf'  DEVICE=PDF;
GOPTIONS FTEXT=ZAPF BORDER ROTATE=LANDSCAPE;
SYMBOL1 V=DOT H=.1 C=BLACK I=JOIN L=1; 
AXIS1 LABEL=('Temps');
AXIS2 LABEL=(a=90 'Serie');
PROC GPLOT DATA=transfo;
  PLOT logZt*NBOBS=1 / haxis=axis1 hminor=0 vaxis=axis2 vminor=0;
TITLE "TRANSFORMATION LOG";
run; 

/* 
De manière plus systématique, on peut considérer la transformation de Box-Cox et utiliser
la macro fournie dans le module SAS/EST  %BOXCOXAR:
*/

%BOXCOXAR(serie, Zt, AR=5, LAMBDALO=-2, LAMBDAHI=2, NLAMBDA=21);

/* 
Remarques:
1.  En principe, à moins de bizarreries extraordinaires, vous devriez décider de:
    i)  Ne pas transformer.
    ii) Prendre l'une des trois transformations populaires: racine, inverse ou log.

2.  La méthode de Box-Cox offre l'avantage d'être plus objective.  Elle repose sur
    une approximation par un AR(p).  Puisqu'on soupçonne que les données sont saisonnières
    et avec une composante saisonnière, il est avantageux de prendre un p assez grand 
    (p=5 et même plus grand, par exemple p=20, tout dépendant de la taille de la série).  
*/

/* =============== Definition d'un nouveau jeu de données =============== */
DATA serie2;
  SET serie;
  Zttrans = Zt;
RUN;

/*
Remarque: Vous remarquerez que par défaut, Zttrans est définie comme la série originale.  Si vous
          décidez de transformer, vous devez l'indiquer dans le jeu de données ci-dessus.  Ceci a
          des implications lors du calcul des prévisions, puisque transformer implique
          un changement d'échelle.  Par exemple, si vous décider de transformer avec logarithmique, vous
          devez prendre l'exponentielle pour revenir dans l'échelle originale.
*/

/* ================================================================= */
/* ======== Étape 3: Retrait d'une composante saisonnière? ========= */
/* ================================================================= */

/*  
La plupart des séries sur CANSIM/CANSIMII affichent une forte composante saisonnière, et il
est habituellement souhaitable de considérer la composante 1-B**12.  Vous devriez considérer
sérieusement cette possibilité.  On considère un nombre de délais élevés afin de regarder ce 
qui se passe autour des délais saisonniers 12, 24, 36,...
*/

PROC ARIMA DATA=serie2;
  IDENTIFY VAR=Zttrans  NLAG=48;
RUN;

/* ================================================================= */
/* ============ Etape 4:  Étude de la stationnarite? =============== */
/* ================================================================= */
/*
L'objectif de cette section est de cerner la question de la stationnarité.
Ainsi, minimalement. est-ce que la série choisie semble stationnaire par rapport à la moyenne?
Si ce n'est pas le cas, il est recommendé de prendre une première différence.  Dans PROC ARIMA,
ceci implique la modélisation de la variable Zttrans(12) = (1-B**12)Zttrans.

L'argumentation se fonde sur un résultat théorique: comme fonction de la taille échantillonnale n, 
c'est-a-dire a mesure que n augmente, pour un délai k fixé, si vous remarquez que les r(k) tendent
dangereusement vers un, vous avez un doute de la stationnarité.  Si vous doutez de la stationnarite, 
vous pouvez considérer de prendre une première différence dite régulière: Zttrans(1) = (1-B)Zttrans.
*/

DATA split1;
  TAILLE = 100;    /* Taille de la série chronologiques */
  SET serie2;
  IF NBOBS < TAILLE/3;
RUN;

DATA split2;
  TAILLE = 100;    /* Taille de la série chronologiques */
  SET serie2;
  IF NBOBS < 2*TAILLE/3;
RUN;

PROC ARIMA DATA=split1;
  IDENTIFY VAR=Zttrans  NLAG=48;
RUN;
PROC ARIMA DATA=split2;
  IDENTIFY VAR=Zttrans  NLAG=48;
RUN;
PROC ARIMA DATA=serie2;
  IDENTIFY VAR=Zttrans  NLAG=48;
RUN;

/*
  Vous pouvez repeter la petite analyse de cette section avec (1-B**12)Zttrans.
*/ 

PROC ARIMA DATA=split1;
  IDENTIFY VAR=Zttrans(12)  NLAG=48;
RUN;
PROC ARIMA DATA=split2;
  IDENTIFY VAR=Zttrans(12)  NLAG=48;
RUN;
PROC ARIMA DATA=serie2;
  IDENTIFY VAR=Zttrans(12)  NLAG=48;
RUN;


/* ================================================================= */
/* ============= Conclusion de l'étape d'identification ============ */
/* ================================================================= */
/*
  À ce stade, vous pensez disposer d'une certaine serie
  stationnaire, et il est temps d'utiliser PROC ARIMA de SAS.  
  Si vous decidez par exemple de considerer (l-B)(1-B^12), alors
  vous devrez en tenir compte dans l'énoncé IDENTIFY de PROC ARIMA.  
  Si une première différence régulière et une première différence saisonnière
  sont requises, alors dans PROC ARIMA vous devrez spécifier Zttrans(1,12).
  Remarque: si deux différences régulières sont requises, vous devez spécifier
            Zttrans(1,1) = (1-B)**2 Zttrans.  Si vous spécifier Zttrans(2), vous obtiendrez
            la différentiation (1-B**2) Zttransf qui n'est pas équivalente à 
            Zttrans(1,1) = (1-B)**2 Zttrans = (1- 2*B + B**2) Zttrans.
*/

ODS PDF CLOSE;