/* exercice 1 */ /* chronique pub (moteur v6) */ /* ventes en fonctions des dépenses de publicité */ title; data pub; set chron6.pub; date=_n_; input=dif(pub); output=dif(ventes); run; proc print data=pub; run; proc sgplot data=pub; series x=date y=ventes; /* en noir */ series x=date y=pub/Y2AXIS; /* en bleu */ run; /* la courbe bleue a l'air de "précéder" la courbe rouge (avec un peu de bonne volonté...), il n'est pas absurde de chercher à expliquer ventes par pub */ proc sgplot data=pub; series x=date y=output; series x=date y=input/Y2AXIS; run; /* on commence par regarder ACF/PACF des deux chroniques pour voir si elles sont stationnaires */ proc arima data=pub; identify var=pub; run; identify var=ventes; run; /* decroissance lente pour pub, on derive */ /* on dérive aussi pour ventes, meme si la chronique a déja l'air stationnaire car on cherche un modele ventes en fonction de pub, et pas ventes cumulées en fonction de pub */ /* comme dans le cours, on cherche à blanchir la chronique (1-B)pub*/ identify var=pub(1); run; estimate p=2 plot; run; /* un filtre grossier suffit pour le préblanchiment (le but est seulement de pouvoir identifier la fonction de transfert) */ /* au passage: q=2 marche mieux que p=2 (cf Portmanteau et ACF/PACF meilleurs) */ /* on peut aussi filtrer par un MA(2), cela ne change rien sur l'interprétation du corrélogramme croisé qui suit */ identify var=ventes(1) crosscor=pub(1); run; /* notez l'avertissement de SAS: les deux chroniques ont été préfiltrées avec AR(2)*/ /* pics aux retards 2 3 4 6 7 (2 et 8 limites)*/ /* CONCLUSION: Intervalle de longueur 6-7 dans les corrélations croisées, et décalage de 2 */ /* vu les oscillations sinusoidales, il doit y avoir un degré au moins 2 au dénominateur */ /* on essaie différentes fonctions de transfert avec dénominateur de degré 2*/ /* numérateur constant: */ estimate input=(2 $ /(1,2) pub) plot; /* 2 avant $: décalage */ run; /* les ACF/PACF des résidus indiquent qu'ils ne sont pas stationnaires (on ne peut pas identifier un ARMA simple: ici longue mémoire) : la fonction de transfert doit etre enrichie */ /* numérateur degré 1: */ estimate input=(2 $ (1)/(1,2) pub) plot; run; /* encore des effets "long-terme" (dans PACF, au décalage de 6)*/ /* on peut donc essayer directement avec un polynome de degré 2 au numérateur*/ /* sinon on essaie à ce stade d'identifier un modèle ARMA sur les résidus */ /* des pics dans le PACF, "décroissance" sinusoide lente de l'ACF */ /* il y a certainement un mix AR+MA */ /* on tatonne... (en verifiant que les coefs sont significatifs) */ /* n'hésitez pas à essayer d'autres modèles et à départager les modèles valides selon SBC, AIC */ /* si on s'arrête à un terme au numérateur, on finit par tomber sur le modèle: */ estimate p=2 q=(2) input=(2 $ (1)/(1,2) pub) plot; run; /* la constante n'étant pas significative: */ estimate p=2 q=(2) input=(2 $ (1)/(1,2) pub) noconstant plot; run; /* Std Error Estimate 4.624289 AIC 561.3748 SBC 579.1778 */ /* avec polynôme de degré 2 au numérateur: */ estimate input=(2 $ (1,2)/(1,2) pub) plot; run; /* il faut identifier le modèle ARMA des résidus (pas évident: il y a toujours des corrélations long-terme) */ estimate p=1 q=1 input=(2 $ (1,2)/(1,2) pub) plot; run; /* pas bon */ /* on passe aux modèles à p+q=3 */ estimate p=2 q=1 input=(2 $ (1,2)/(1,2) pub) plot; run; /* presque bon pour ACF et PACF, mais Portmanteau pas bon */ estimate p=1 q=2 input=(2 $ (1,2)/(1,2) pub) plot; run; /* pas bon */ /* modèles à p+q=4 */ estimate p=1 q=3 input=(2 $ (1,2)/(1,2) pub) plot; run; /* pas bon */ estimate p=3 q=1 input=(2 $ (1,2)/(1,2) pub) plot; run; /* pas bon */ estimate p=2 q=2 input=(2 $ (1,2)/(1,2) pub) plot; run; /* ok */ /* s=4.35 AIC=552.6 SBC=578 */ /* mais le coef de décalage 1 du MA et la constantes ne sont pas significatifs */ /* le modèle suivant est légèrement meilleur */ estimate p=2 q=(2) input=(2 $ (1,2)/(1,2) pub) plot; run; /* s=4.33 AIC=550 SBC=573 */ /* tout est ok, normalité des résidus correcte */ /* on garde le meilleur (celui-ci)*/ /* modèle: (1-B)Ventes_t = -0.137 + (0.879 + 0.828B - 0.573B^2) / (1 - 1.251B + 0.745B^2) (1-B)Pub_{t-2} + (1 - 0.60891B^2) / (1 - 1.117B + 0.757B^2) eps_t */ /* REMARQUE: la constante n'étant pas significative, mais si on l'enlève: */ estimate p=2 q=(2) input=(2 $ (1,2)/(1,2) pub) noconstant plot; run; /* s=4.306363 AIC 548.897 SBC 569.2433 */ /* ATTENTION: sur certaines version de SAS, il y a instabilité nummérique */ /* modèle souvent proposé dans la littérature, dont on peut se satisfaire: estimate p=2 q=2 input=(2 $ (1,2)/(1,2) pub) plot; run; */ forecast lead=12 back=12 id=date out=b; run; /* même remarque que précédemment: pour une "vraie" prévision, il faudrait des données "neuves" pour pub */ proc sgplot data=b; band Upper=u95 Lower=l95 x=date / LegendLabel="95% Confidence Limits"; scatter x=date y=ventes; series x=date y=forecast; run; /* exercice 2 */ /* Hudson */ /* il s'agit du nombre de peaux de rats musqués et visons vendues par la compagnie de la Baie d'Hudson entre 1848 et 1909 */ /* rats musqués: proies visons: prédateurs on voudrait expliquer le nombre de rats musqués en fonction du nombre de visons */ proc print data=chron6.hudson; run; data hudson; set chron6.hudson; date=1848+_n_-1; lvisons=log(visons); lmusques=log(musques); run; proc sgplot data=hudson; series y=visons x=date; series y=musques x=date/y2axis; run; /* variance pas stable: passage au log */ proc sgplot data=hudson; series y=lvisons x=date; series y=lmusques x=date/y2axis; run; /* séries avec peu de données, donc intervalles de confiance assez larges: c'est difficile de juger de la significativité */ proc arima data=hudson; identify var=lmusques nlag=25; run; identify var=lvisons nlag=25; run; /* la periodicité n'entraine pas de pic significatif, donc on choisit de ne pas faire de differentiation saisonniere Par contre, on dérive à l'ordre 1 car l'aspect des graphes n'est pas très stationnaires */ proc arima data=hudson; identify var=lmusques(1) nlag=25; run; identify var=lvisons(1) nlag=25; run; estimate p=(10) plot; run; /* ( pour pré-blanchir visons(1) ) */ identify var=lmusques(1) crosscor=lvisons(1); run; /* les corrélations croisées pour h<0 sont limites, les corrélations croisées pour h=0 et h=1 sont significatives */ /* même démarche que plus haut */ estimate input=(0 $ (1) / lvisons) plot; run; estimate input=(0 $ (1) / lvisons) noconstant plot; run; estimate input=(0 $ (1 2) / lvisons) plot; run; /* les résidus forment un b.b.g. */ /* Std Error Estimate 0.256249 AIC 10.62333 SBC 18.93348 */ /* meilleur */ /* avec un dénominateur: */ estimate input=(0 $ /(1) lvisons) plot; run; /* difficile de mettre un ARMA sur les résidus (mémoire longue), on met aussi un terme sur le numérateur */ estimate input=(0 $ (1)/(1) lvisons) plot noconstant; /* constante non significative */ run; /* Std Error Estimate 0.247408 AIC 5.588993 SBC 11.87203 */ /* meilleur */ forecast lead=10 back=10 id=date out=b; run; data b; set b; visons=exp(lvisons); musques=exp(lmusques); eforecast=exp(forecast+std*std/2); el95=exp(l95); eu95=exp(u95); run; proc gplot data=b; plot residual*date; symbol i=needle width=1; run; proc sgplot data=b; band Upper=eu95 Lower=el95 x=date / LegendLabel="95% Confidence Limits"; scatter x=date y=musques; series x=date y=eforecast; run; /* ne pas oublier qu'il y a finalement peu d'observations pour construire le modèle */