SARIMA ja kampanjan vaikutuksen arvioiminen R-ohjelmalla

Käsittelen nyt SARIMA:a ja kampanjan vaikutuksen arvioimista. Haluan tietää, onko mainoskampanjalla ollut vaikutusta myyntiin. Käyttämässäni harjoitusaineistossa kampanja on pidetty aikavälillä 01/01/2009- 31/03/2009. On hyvä muistaa, että mainoskampanjan vaikutus voi ulottua kampanjan kestoa pidemmälle aikavälille. Esimerkiksi ihmiset saattavat ostaa tuotetta vasta kuukauden päästä kampanjan loppumisesta, kun mainos palautuu mieleen. Ajatellaan mainoskampanja tässä vain mainostamiseksi, eli oletetaan, ettei meillä ole kampanjahintaa. Tällä tavoin voimme mitata myyntiä ja olettaa myös, että kohonneet myyntiluvut tarkoittavat uutta tuloa, koska myyntikatteet pysyvät samoina kuin ennen kampanjaa. Tehtävänä on nyt saada selville, onko kampanja ollut hyödyllinen, eli onko se lisännyt myyntiä (kampanjan tarkempi hyöty selviää tietysti vasta, kun sen kokonaiskustannukset ovat tiedossa) Alla olevasta kuvasta näemme yrityksen myynnit sekä kampanjan aikavälin. Pyrin ensin ennustamaan myynnit kampanjan aikavälille, jos kampanjaa ei olisi pidetty. Eli muodostan SARIMA-mallin ja ennustan myynnit perustuen puhtaasti myyntihistoriaan. Tämän jälkeen vähennän ennustetut myynnit toteutuneista ja näemme, onko kampanja ollut hyödyllinen.

Käytän ennustamisessa Rob J Hyndman:n suosittelemaa tapaa, jossa jaan datan ns. Training periodiin ja Test periodiin. Training periodissa muodostan mallin ja ennustan sillä Test periodin arvoja vastaan. Testaan siis mallia oikeaan dataan ennen kuin ennustan sillä halutulle aikavälille. Tällä tavoin en ennusta täysin tuntemattomaan, vaan mallini on onnistunut ennustamisessa jo aiemmin, eli on ennustanut parhaiten Test periodin aikana. Hyndman neuvoo, että Test periodiksi tulisi valita vähintään yhtä pitkä aikaväli kuin ennustusaikaväli on, eli jos ennustamme 12 kuukautta, niin Test periodin on oltava vähintään 12 kuukautta. On tärkeää huomata, että ns. Forecast accuracy:a, eli ennustustarkkuutta mitataan tällä tavoin eikä sillä, miten hyvin malli sopii Training periodin dataan. Tosin voimme mallia valitessamme hyödyntää malleja, joissa on alhaisin Akaike Information Criteria (AIC)-arvo, eli ikään kuin esikarsimme malleja AIC-arvon perusteella. AIC-arvo kertoo mallien suhteellisesta paremmuudesta, eli malli, jossa AIC on pienin on suhteellisesti paras malli, koska informaatiokadon estimaatti on pienin. Hyndman:n mukaan AIC-arvoa ei voi käyttää mallien vertailuun, jos differensointien määrä vaihtelee, eli AIC-arvoa voi vertailla vain jos differensoinnit ovat samat. Esimerkiksi voimme vertailla malleja, joissa on yksi kausidifferensointi keskenään.

Dataa katsomalla huomaamme, että meidän on kokeiltava ns. Box-Cox-transformaatiota, eli muunnamme aineiston logaritmi-muotoon. Tällä tavoin varianssi pysyy aikavälillä suhteellisen vakiona, joka on tärkeä oletus ARIMA:a/SARIMA:a, eli ns. Additiivisia malleja käytettäessä. Alla olevassa kuvassa näkyvät myynnit logaritmi-muunnoksen jälkeen. Vaikka alkuperäisessä muodossa myynnit kasvavat absoluuttisesti enemmän aikasarjan-lopussa, niin silti suhteellisesti kasvu on melko vakaata, joten logaritmi-muunnos on hyödyllinen.

Seuraavaksi voimme katsoa ACF (autokorrelaatio) ja PACF (osittaisautokorrelaatio) kuvat, jotta näemme, onko aikasarjassa yksikköjuurta (Unit root). Jo kuvia katsomatta on selvää, että sellainen on, koska keskiarvo vaihtelee ajanjakso aikana, eli esim. 1999-2000 välillä mitattu keskiarvo ei ole sama kuin vuosien 2004-2005 välillä mitattu. Meidän pitää muokata aineistoa differensoimalla, jotta pääsemme yksikköjuuresta eroon. Kuitenkin ennen differensointia on huomattava, että ainestossa on selkeä nouseva trendi, ja jos haluamme trendin mukaan malliimme, on meidän huomioitava se ennen differensointia, sillä differensointi poistaa trendin (Ajattele esim. kasvavaa trendillistä aikasarjaa 5,10,15,20, josta otamme peräkkäisten arvojen differensoinnin. Saamme sarjaksi 5,5,5 (10-5,15-10,20-15) eli aikasarja muuttuu vaakaviivaksi) . Voimme ottaa kausidifferensoinnin (Seasonal differencing) ja/tai differensoinnin (Differencing) peräkkäisille arvoille. Kausidifferensointi tarkoittaa, että vähennämme jokaisesta arvosta sen aiemman vuoden arvon: Esim. joulukuun 2001 kausidifferensoitu arvo on joulukuun 2001 arvo "miinus" joulukuun 2000 arvo. Peräkkäisten arvojen differensointi taas tarkoittaa, että joulukuun 2001 arvo saadaan vähentämällä siitä marraskuun 2001 arvo.

Jaan nyt aineiston "Trainingset":iin 01/1999-12/2007) ja "Testset":iin (01/2008-12/2008), jotta voin arvioida ennustustarkkuutta, josta edellä puhuttiin. Alla olevassa kuvassa ovat trainingsetin ACF- ja PACF -kuvat

ACF-kuvasta näemme vahvaa autokorrelaatiota monelle viivelle. Alla olevassa kuvassa ovat kuukausittaiset keskiarvot Trainingsetille ja näemme, että keskiarvot vaihtelevat kuukaudesta riippuen. Heinä-ja elokuu ovat selkeästi myynnillisesti parempia kuukausia kuin muut. (Kun kausivaihtelua ei esiinny, niin kuvan keskiarvot olisivat likimain samalla tasolla. )

Teen ensimmäiseksi kausidifferensoinnin ja katson sen jälkeen, onko toinen differensointi tarpeen. Kausidifferensointi kannattaa tehdä aina ensin, koska voi olla, että jo se riittää, eli toista peräkkäisten arvojen differensointia ei tarvita (Differensoinnissa katoaa informaatio (Trendi)). Tavoite on saada aineisto stationaariseen muotoon, jotta mallista saadut tulokset olisivat luotettavia. Ei-stationaarinen aineisto altistaa mallin tulokset autokorrelaatiolle, joka vääristää keskivirheet, eli esim. luottamusvälit eivät olisi luotettavia. Stationaarisessa aikasarjassa keskiarvo ja varianssi eivät vaihtele ajasta riippuen, joten sitä on mahdollista ennustaa. Stationaarisessa aikasarjassa muuttujan regressiokerroin on vakio ajasta riippumatta, kun taas yksikköjuurellisessa aikasarjassa regressiokertoimen tulisi muuttua ajassa, jotta se kuvaisi aikasarjaa oikein. Ei-stationaarisessa aikasarjassa (yksikköjuurellisessa) voisi olla esim. pitkiä nousuja ja tämän jälkeisiä pitkiä laskuja. Kuinka voisimme ennustaa tai kuvata tällaista aikasarjaa yhdellä regressiokertoimella? Stationaarisuus-ehto on tärkeä myös, koska silloin Law of the large numbers- ja Central limit theorem -lait pätevät.

Alla olevassa kuvassa on kausidifferensoitu aikasarja, joka ei ADF-testin mukaan sisällä yksikköjuurta. KPSS-testi osoittaa sarjan olevan nyt stationaarinen, vaikka asia ei silmämääräisesti aivan selvä olekaan.

Nyt kun aikasarja on stationaarisessa muodossa, voin katsoa ACF-ja PACF-kuvat, jotta saan informaatiota mallin valintaan. Yleistäen voisi sanoa, että ACF-kuvasta katsotaan tarvittavat MA-termit ja PACF-kuvasta tarvittavat AR-termit. Esimerkiksi ACF-kuvassa ensimmäisen viiveen autokorrelaatio viittaisi siihen, että MA-termi sopisi malliin. ACF-kuvan 12. viiveen kohdalla oleva piikki voi tarkoittaa, että MA-termi tulee lisätä myös kausivaihtelukomponenttiin. Lisää tietoa mallinvalinnasta löytyy Robert Nau:n sivuilta: http://people.duke.edu/~rnau/seasarim.htm

Voin myös käyttää Auto.arima-toimintoa Forecast-pakkauksesta. Auto.arima kertoo oman ehdotuksensa parhaaksi malliksi. Auto.arima ehdottaa kausidifferensointia+trendiä käytettäväksi. Auto.ariman ehdotus on tarkalleen ARIMA(2,0,0)(2,1,0)[12]+Drift. Ensimmäinen luku vasemmalla on AR-termien määrä, keskimmäinen on peräkkäisten arvojen differensointien lukumäärä, kolmas luku taas on MA-termien määrä. Näiden vieressä ovat vastaavat arvot kausivaihtelukomponentille. [12] kuvaa kausivaihtelun periodien lukumäärää.  Jos data olisi neljännesvuosittain julkaistua, niin periodien lukumäärä olisi 4.

Näytän seuraavaksi, miten trendin (Drift) lisääminen vaikuttaa mallin ennusteisiin. Alla olevassa kuvassa ovat ARIMA(2,1,0)(2,1,0)[12] ja ARIMA(2,0,0)(2,1,0[12]+Drift. Driftin lisäys on sama asia kuin lisäisi aikamuuttujan malliin, eli muuttujan joka alkaa arvosta 1 ja kasvaa yhdellä numerolla aina sarjan loppuun saakka. Esim 12 kuukauden pituisen aikasarjan aikamuuttuja voisi olla numerot 1-12. Trendin lisäys tarkoittaa sitä, että mallin ennusteet ovat trendin mukaisia. Se onko tämä hyvä vai huono asia on itse päätettävä. Hyndmann varoittaa trendin käytöstä, koska luottamusvälit ovat kapeita, kuten alla olevasta kuvasta näemme. Trendin käyttö mallissa voi tosiaan tuottaa liian optimistisen kuvan tulevaisuuden myynneistä.

Malli ilman trendiä (yläpuolella) ja malli trendin kanssa (alapuolella)

Nyt, kun aikasarja on stationaarisessa muodossa ja olen nähnyt ACF-ja PACF-kuvat sekä Auto.ariman mallisuosituksen voin valita malleja ja testata niitä testidataan, eli "Testset":n, jonka alussa tein. Ennustan siis eri malleilla aitoon vuoden 2008 dataan ja vertaan mallien tarkkuutta.
Valitsen tarkasteluun seuraavat mallit:
Malli1: ARIMA(2,0,0)(2,1,0)[12]+Drift
Malli2: ARIMA(2,0,0)(2,1,1)[12]+Drift
Malli3: ARIMA(2,0,0)(1,1,1)[12]+Drift
Malli4: ARIMA(2,0,2)(2,1,0)[12]+Drift
Malli5: ARIMA(1,0,2)(2,1,0)[12]+Drift

Seuraavaksi testaan, kuinka tarkasti mallit ennustavat vuoden 2008 arvot. R antaa erilaisia arvoja, kuten ME, RMSE, MAE,MPE,MAPE ja MASE. Malleista tarkoin on se, jolla kyseiset arvot ovat pienimmät. Ei ole kuitenkaan varmaa, että löytyisi malli, jossa kaikki arvot olisivat toisia pienempiä. Alla kuvassa mallit 1 ja 4 sekä niiden ennusteet vuodelle 2008

Kuten kuvista näkee on ennusteiden tarkkuudessa toivomisen varaa, eli kumpikin malli yliestimoi tulevia arvoja. Tätä selittää trendin mukana olo, eli kokeilen seuraavaksi samoja malleja 1 ja 4, kun trendi on otettu pois. Huom. Trendin/driftin pois ottaminen ei tarkoita, että ennusteet eivät edelleen kasvaisi, sillä kausivaihtelukomponentti sisältää AR-termin, eli SAR-termin (SAR=Seasonal AR). Jos haluat tietää enemmän SAR- ja SMA -termien vaikutuksesta malliin, niin katso sivun alaosan kappale SAR ja SMA.

Testaan vielä nämä mallit, kun olen lisännyt niihin SMA-termin (Seasonal MA-termi). SMA-termi voisi sopia malliin, koska ACF-kuvassa oli negatiivinen piikki 12-lagin kohdalla (ks. Rule 13: http://people.duke.edu/~rnau/arimrule.htm).

Nyt ennusteet osuvat jo paremmin. Parhaiten onnistuu viimeinen malli: ARIMA(2,0,2)(2,1,1)[12]. Varmistan vielä, että tällä mallilla ei esiinny residuaalien autokorrelaatiota. Käytän testaamiseen Ljung-Box-testiä ja testi ei hylkää nollahypoteesia, joten virhetermeissä ei esiinny autokorrelaatiota. Käytän siis tätä mallia.

Seuraavaksi ennustan valitulla mallilla ARIMA(2,0,2)(2,1,1)[12] kaksi vuotta eteenpäin alkaen 01/2009, eli saamme ennusteet kampanjan ajalle ja kampanjan jälkeiselle ajalle. Olen kiinnostunut myös kampanjan jälkeisestä ajasta, koska kampanjan vaikutus on voinut kestää pidempään kuin kampanja.

Kuvasta näemme, että koko kampanjan ajan (01/2009-03/2009) toteutuneet myynnit ovat ennustettuja arvoja korkeammat. Tämä kertoo siitä, että kampanja on onnistunut ainakin myynnillisesti (Ottamatta kantaa kampanjan kokonaiskustannuksiin). Mielenkiintoinen seikka on, että kampanjan vaikutus näyttää ulottuvan vuoden 2009 loppuun saakka. Päättymiskohdan voi katsoa kohdasta, jossa luottamusväli leikkaa toteutuneet myynnit. Alla toteutuneiden myyntien ja ennustettujen myyntien (point forecasts) erotukset vuodelle 2009.

Jan   Feb   Mar   Apr   May  Jun   Jul   Aug   Sep   Oct   Nov   Dec
61.5 61.0  74.2  74.5   88.4   69.2 92.2 99.6   86.0  83.4  86.6  95.6

Alla on vielä toteutuneiden myyntien ja ennustettujen myyntien 95%:n luottamusvälin ylärajan erotus vuodelta 2009. Luottamusväli sisältää ylä-ja alarajan, joiden välissä piste-ennustus on.Huom. Piste-ennusteen luku ei sinällään ole yhtään luotettavampi luku kuin muukaan yksittäinen luku luottamusväliltä. Luottamusväli kertoo meille, että 95% kerroista, kun teemme luottamusvälin saamme oikean arvon luottamusvälillemme. Oikea luku voi olla siis mikä tahansa arvo luottamusvälillä, niinpä on varminta laskea kampanjan vaikutus luottamusvälin ylärajalta.

 Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
35.2  30.8  34.9  31.7   40.8    8.1  19.7  22.0   17.8  20.8  29.8  30.3




SAR- ja SMA

SAR-termin vaikutus selviää parhaiten, kun ajattelemme kahden mallin eroa, joista toinen on ARIMA(0,0,0)(0,1,0)[12]-malli ja toinen on malli, jossa on lisäksi yksi SAR-termi  Malli on siis: ARIMA(0,0,0)(1,1,0)[12]
Estimoin nyt molemmat mallit Trainingperiodin dataan.

Ensimmäisen mallin ennuste tulevalle kuukaudelle on aina vuoden takaisen kuukauden arvo, eli tammikuun 2008 ennustettu arvo on tammikuun 2007 arvo. Mallissa ei ole siis nousevaa tai laskevaa trendiä, vaan sen ennusteet ovat täysin kopioita ennustetta edeltävältä 12 kuukaudelta.
Kun estimoimme toisen mallin Trainingperiodin dataan, saamme SAR(1) arvon, joka on 0.7237. Nyt tammikuun 2008 ennuste on: tammikuun 2007 arvo "plus" 0.7237*(tammikuun 2007 arvo "miinus" tammikuun 2006 arvo). Eli koska datassa on nouseva trendi (tammikuun arvo 2007 > tammikuun arvo 2006 jne.), jonka takia uudemman kuukauden arvot ovat aina vuoden takaista korkeammat, niin myös mallin ARIMA(0,0,0)(1,1,0)[12] ennusteissa on nouseva trendi, joskin se on ajassa loiveneva (SAR(1)-kerroin 0.7237 < 1 aiheuttaa sen, että ennusteiden välinen ero kapenee jatkuvasti, kunnes eroa ei lopulta enää ole ja ennusteet pysyvät vakiona). SAR-termin johdosta malliin saadaan siis mahdollinen trendi mukaan. On huomattava, että mallissa "trendi" ei ole sama kuin deterministinen aikatrendi (mallit 1-5 joissa drift mukana) koko ajanjaksolle, vaan se on riippuvainen edellisten kahden vuoden takaisten kuukausien erotuksesta, eli esimerkissämme tammikuun 2008 arvo riippuu tammikuun 2007 ja tammikuun 2006 arvojen erotuksesta. Jos lisäisimme SAR(2)-termin malliin, niin tammikuun 2008 arvo riippuisi lisäksi myös tammikuun 2006 ja tammikuun 2005 arvoista. Voisi ajatella, että SAR(3) ja suuremmat SAR(x):t ovat sen verran harvinaisia, että niitä tuskin malleihin koskaan tarvitsee. SAR(3) esim. aiheuttaisi sen, että mallin ennuste tammikuulle 2008 olisi riippuvainen tammikuun 2005 ja tammikuun 2004 arvoista. Alla vielä kuvat kahdesta mallista.

SMA-termin vaikutusta voidaan tarkastella mallilla ARIMA(0,0,0)(0,1,1)[12]. Estimoimalla mallin Trainingperiodin dataan saadaan SMA(1)-kerroin, jonka arvo on 0.472. Nyt voimme laskea ennusteen tammikuulle 2008 seuraavasti: tammikuun 2007 arvo "plus" 0.472*(tammikuun 2007 oikea arvo "miinus" tammikuun 2007 ennustettu arvo). Jos siis ennuste on mennyt alakanttiin, niin ennuste tammikuulle 2008 on suurempi kuin tammikuun  2007 arvo. Tässä mallissa SMA(1)-termi vaikuttaa vain ensimmäisen vuoden ennusteisiin, koska SMA-termi tarvitsee vaikuttaakseen myös ennustetun arvon. Toisaalta vaikutus on tavallaan pidempi, jos mukana on AR tai SAR-termejä, koska nämä luovat uusia arvoja vanhoista arvoista, joihin SMA on vaikuttanut, eli ennusteet esim. mallien ARIMA(1,0,0)(0,1,0)[12] ja ARIMA(1,0,0)(0,1,1) välillä eroavat myös, kun SMA(1)-termin vaikutus on loppunut. Alla kuva ARIMA(0,0,0)(0,1,1)[12]-mallista.