Logistinen regressio ja luottoriskien estimointi R-ohjelmalla

Tässä tekstissä käsittelen logistista regressiota ja luottoriskien estimointia. Muodostan mallin, jolla pankin asiakkaiden maksukyvyttömyyden (default) todennäköisyyttä voidaan arvioida perustuen erilaisiin muuttujiin. Muuttujina ovat mm. asiakkaan luottohistoria, haetun luoton määrä, työllisyystilanne, säästöt, asumismuoto ja ikä. Data (german.data) ja lisätietoa datasta löytyy täältä. Menetelmät, joita käytän perustuvat D, Sharman oppaaseen, joka löytyy täältä. Datasta Default-muuttuja toimii selitettävänä muuttujana ja muut 20 muuttujaa selittävinä muuttujina.
Default-muuttuja on binäärimuuttuja, joka saa arvon 1 jos asiakas on maksukyvytön ja arvon 0 jos asiakas ei ole maksukyvytön, eli pystyi maksamaan luottonsa takaisin. Tarvitsemme siis luotonmyöntöön vaikuttavien tietojen lisäksi tiedon siitä, kuinka asiakkaille kävi, eli pystyivätkö he maksamaan luottonsa vai eivät. Tämän tiedon avulla voimme muodostaa mallin ja mahdollisesti tehdä havaintoja, kuten esim. että työssäkäyvät selviytyvät luotoistaan työttömiä useammin. Jos nyt kaksi muuten täysin identtistä asiakasta tulisi hakemaan pankista lainaa, niin mallimme mukaan työllinen selviytyisi luotostaan todennäköisemmin. Se, että toinen maksaa lainansa takaisin todennäköisemmin ei tarkoita välttämättä, että toinen ei saisi lainaa ollenkaan. Lainanantaja voi nimittäin korkoa kasvattamalla pitää tuotto-odotuksensa vakiona, eli vaikka maksukyvyttömyyttä esiintyisi enemmän tietyllä ryhmällä, niin suuremmat korot voivat kompensoida nämä tappiot. Alla on esimerkki yksinkertaisesta yhtälöstä, joka kuvaa riskillisen tuoton ja riskittömän tuoton yhteyttä. Riskitön tuotto "r(g)" voi kuvata esim. EKP:n talletuskorkoa, eli tuottoa, jonka pankki saa rahalleen ilman riskiä. Yhtälössä "(1-p)" on "ei maksukyvyttömyyden" todennäköisyys, "u" on raha, joka saadaan takaisin maksukyvyttömyyden sattuessa (oletetaan, että kaikki menetetään, jos asiakas on maksukyvytön, eli u=0), "z" riskipreemio ja "r(f)" riskillinen tuotto . Riskipreemio kuvaa tuottoa, jonka pankki haluaa lisäksi, kun riskillisen tuoton odotusarvo vastaa riskitöntä tuottoa. Yhtälöstä (1) voidaan ratkaista ns. spreadi, joka kuvaa riskillisen- ja riskittömän tuoton erotusta (S=r(f)-r(g)). Tämä spreadi on esitetty yhtälössä (2). (Huom. ratkaisussa on oletettu, että p ja r(f) ovat pieniä, jotta p*r(f)~0

1.  (1-p)(1+r(f)) + p*u = 1+r(g)+z
2.  S=p+z

Spreadiin vaikuttaa siis sekä maksukyvyttömyyden todennäköisyys että riskipreemio. Vielä yksinkertaistan, kun oletetaan, että riskipreemio on nolla, voimme kirjoittaa spreadin S=p eli r(f)-r(g)=p

Jos siis pankki ei pyydä riskipreemiota ja haluaa odotusarvoisesti saada saman tuoton sekä työlliseltä että työttömältä, niin sen pitää tietää p mahdollisimman tarkasti. Jos työllisen todennäköisyys maksukyvyttömyyteen olisi 1%:n ja työttömän vastaava luku 3%:a, niin pankki pyytäisi esimerkin tilanteessa työlliseltä korkoa ¨~1%+riskitön korko ja vastaavasti työttömältä  ~3%+riskitön korko. Tällä tavoin pankin odotusarvo kummallekin saatavalla olisi sama.


Nyt takaisin malliin ja sen luomiseen. Ensin lataan R:ään datan ja nimeän sarakkeet, koska niillä ei ole valmiina nimiä.
R-koodi:

rm(lists=ls())    #poistaa aiemmat tallennetut muuttujat. Hyvä käyttää, jos attach-komentoa käytetään.
data <- read.delim(file.choose(), header=FALSE,sep="")

names(data) <- c("check_acc", "duration", "history", "purpose", "amount", "sav_acc", "employment", "install_rate", "ptatus", "other_debtor", "resid", "property", "age", "other_install", "housing", "other_credits", "job", "num_depend", "teleph", "foreign", "default")
attach(data)
str(data)                              #str-komennolla nähdään ovatko muuttujat kategoriallisia vai jatkuvia

Seuraavaksi voimme katsoa dataa kuvien avulla. Jatkuvia muuttujia, kuten age voidaan katsoa esim. plot-ja hist-komennoilla. Muuttujia, jotka ovat faktoreina (esim.history) voidaan katsoa lattice-pakkauksen histogram-komennolla. Katsotaan seuraavaksi age-muuttujaa:
R-koodi:

library(lattice)

histogram(age,type="percent",breaks=seq(10,80,by=10))   #breaks-komennolla jaetaan eri ikäryhmiin

Näemme, että aineistossa suurin osa on 20-40-vuotiaita. Voimme tarkastella eri ikäluokkien osuuksia default-luokissa vcd-pakkauksen mosaic-komennolla. Jaetaan age-muuttuja ensin kolmeen ikäluokkaan ja muutetaan default-muuttuja faktoriksi, joka saa arvot 0 ja 1. Sen arvo 1 tarkoittaa maksukyvyttömyyttä. Sitten voimmekin tehdä Mosaic-kuvan.
R-koodi:

agec <- as.factor(ifelse(age <= 30, "0-30",ifelse(age <=40,"31-40","41-80")))
data$default <- as.factor(ifelse(data$default==1,"0","1")) #päivitetään myös data.frameen
attach(data)
library(vcd)
mosaic(~default+age,highlighting="age",main="age & default")

Kuvasta nähdään, että 0-30-vuotiaiden osuus on maksukyvyttömistä suurin. Seuraavaksi tutkin, mitkä ovat maksukyvyttömien osuudet eri ikäluokissa. Tämän voi tehdä table-komennolla seuraavasti: 
R-koodi:

defaultagetable <- table(default,age) # Defaultin luokat menevät table:n riveiksi ja age:n luokat sarakkeiksi
prop.table(defaultagetable,margin=2) # Lasketaan maksukyvyttömien ja maksukykyisten osuudet
a <- barplot(prop.table(defaultagetable,margin=2),beside=T,legend.text=c("Maksukykyisten osuus",
"Maksukyvyttömien osuus"),yaxt="n",ylab="%",xlab="age") # Tehdään pylväskuvio
axis(2,at=seq(0,1,by=0.1))
text(a, 0, round(fdf,2),cex=1.1,pos=3)

 Testaan nyt Khiin neliö -testillä (Chi-square test), ovatko nämä erot ikäluokkien välillä tilastollisesti merkitseviä. "0-30" -luokka eroaa muistaa luokista merkitsevästi, kun taas "31-40" ja "41-80" -luokat eivät eroa keskenään merkitsevästi.

Seuraavaksi lisään "age & default" -kuvaan kolmannen muuttujan, joka on other_debtor, eli muuttuja, joka on kuvattu aineistokuvauksessa seuraavasti:

Other debtors / guarantors 
A101 : none 
A102 : co-applicant 
A103 : guarantor 


Muuttuja siis kertoo, onko hakijalla lainalle takaaja ja hakeeko hakija lainaa yksin. Katsotaan seuraavaksi kuvaa.
R-koodi:

mosaic(default~age+other_debtor,main="age, default & other_debtor")

Kuvasta nähdään, että jokaisessa ikäryhmässä "co-applicant"-hakijat olivat suhteellisesti useimmin maksukyvyttömiä. Hakijat, joilla oli takaaja, pystyivät joka ikäluokassa suhteellisesti useimmin maksamaan lainansa takaisin. 

Seuraavaksi muutan osan muuttujista faktori-muotoon, koska ne eivät ole numeerisia muuttujia eli ne saavat arvoja esim. välillä 1-4 ja kuvaavat luokittelua. 
R-koodi:

data$install_rate <- as.factor(data$install_rate)
data$resid <- as.factor(data$resid)
data$other_credits <- as.factor(data$other_credits)

data$num_depend <- as.factor(data$num_depend)

Nyt  jaan datan "trainingdataan" ja "testdataan", eli jaan datan kahteen osaan, joista toisesta (trainingdata) muodostan mallin, jolla ennustan toisen arvoja (testdata). Tosiasiassa mallin toimivuus tulisi testata myös cross-validationilla ja/tai jakamalla data kolmee osaan, joista kolmatta validationsettiä ei käytetä mallia luodessa, vaan sitä kokeillaan vain kerran, kun malli on valmis. Kun mallia verrataan testdataan, niin muuttujia lisäillessä ja mallia "tuunatessa" voi syntyä ns. overfitting, jolloin malli toimii testdataan paremmin kuin mallin ulkopuoliseen validationsettiin. Tavoitteena on siis muodostaa malli, johon syöttämällä testdatan muuttujien arvot voidaan saada mahdollisimman hyviä ennusteita maksukykyisyydestä. Testdatastahan voimme tarkistaa, kuinka hyvin malli on onnistunut, koska testdatassa on default-muuttuja, joka kertoo maksoiko asiakas lainansa takaisin vai ei (tosiasiassa 1-arvo voi tarkoittaa myös, että asiakas on maksanut lainansa myöhässä tai lainanmaksussa on ollut muita vaikeuksia). Otan trainingdataan sattumanvaraisesti 60%:a arvoista, eli testdataan jää 40%:a arvoista. Huom. Datassa on tällä hetkellä mukana myös jatkuvia muuttujia, kuten age ja amount, joten teen aluksi mallin normaalilla datalla. Kokeilen tämän jälkeen mallia, jossa on pelkästään kategoriallisia muuttujia, eli myös muuttujat, kuten age ja amount on muutettu kategorioihin. Sharma suosittelee käytettäväksi mallissa vain kategorisia muuttujia, mutta toisaalta informaatiokadon vuoksi myös jatkuvien muuttujien malleilla on kannattajansa. Jatkuvan muuttujan muuttaminen kategorialliseksi voi olla hyödyllistä, jos muuttujan yhteys selitettävään muuttujaan ei ole lineaarinen. Epälineaarisia yhteyksiä voi toisaalta mallintaa myös esim. neliöjuuri-muunnoksella, jolloin mallissa käytetään muuttujaa sekä normaalissa että neliöjuuri-muodossa. Epälineaarisuuksia tarkastellessa on aina tarkasteltava muuttujaa myös eri kategorioissa, koska epälineaarisuus saattaa johtua  muuttujan erilaisesta käyttäytymisestä jossain kategoriassa (yhteyden ollessa lineaarinen muissa kategorioissa). Jos muuttuja käyttäytyy eri tavoin eri kategorioissa, on interarction-termien käyttö suositeltavaa. Suosittelen, että muuttujia tarkastellaan aina ennen mallin tekemistä huolella ja mahdollinen kategorisointi on pystyttävä perustelemaan teorialla. Lähtökohtaisesti pyrkisin pitämään jatkuvan muuttujan jatkuvana. Joskus selitettävän ja selittävän -muuttujan välillä löytyy yhteys vasta esim. logaritmi-muunnoksen jälkeen, joten eri muunnokset on hyvä käydä läpi ennen kuin lähtee kategorisoimaan muuttujaa.
R-koodi:

e <- sample(nrow(data), 0.6*nrow(data))  #ottaa sattumanvaraisesti 60%:a rivinumeroista
trainingdata <- data[e,]                              # valitaan saadut rivinumerot trainingdataan ja testdataan

testdata <- data[-e,]

Nyt trainingdatasta voidaan muodostaa logistinen regressiomalli. 
R-koodi: 

model1 <- glm(default~check_acc+duration+history+purpose+
amount+
sav_acc+
employment+
install_rate+
ptatus+
other_debtor+
resid+
property+
age+
other_install+
housing+
other_credits+
job+
num_depend+
teleph+

foreign,family=binomial(link=logit),data=trainingdata)
summary(model1)

Mallissa 1 on useampia tilastollisesti merkitseviä muuttujia. Alla tuloksia. (Standard error suluissa)

Results of the model1

	Dependent variable:
	default
check_accA12	-0.550*
	(0.282)
check_accA13	-0.933**
	(0.461)
check_accA14	-1.944***
	(0.321)
duration	0.036***
	(0.012)
historyA31	0.849
	(0.727)
historyA32	-0.059
	(0.565)
historyA33	-0.240
	(0.613)
historyA34	-0.527
	(0.557)
purposeA41	-1.859***
	(0.496)
purposeA410	-2.696**
	(1.134)
purposeA42	-0.553
	(0.346)
purposeA43	-0.809**
	(0.319)
purposeA44	-2.336
	(1.480)
purposeA45	-0.166
	(0.650)
purposeA46	0.092
	(0.563)
purposeA48	-0.696
	(1.465)
purposeA49	-0.297
	(0.423)
amount	0.0001
	(0.0001)
sav_accA62	-0.541
	(0.377)
sav_accA63	-0.532
	(0.502)
sav_accA64	-1.364*
	(0.734)
sav_accA65	-1.044***
	(0.365)
employmentA72	0.759
	(0.637)
employmentA73	0.454
	(0.617)
employmentA74	0.070
	(0.654)
employmentA75	0.517
	(0.614)
install_rate2	0.150
	(0.389)
install_rate3	0.332
	(0.439)
install_rate4	0.731*
	(0.387)
ptatusA92	-0.084
	(0.490)
ptatusA93	-0.630
	(0.477)
ptatusA94	0.296
	(0.579)
other_debtorA102	0.331
	(0.637)
other_debtorA103	-0.844*
	(0.505)
resid2	1.186***
	(0.400)
resid3	0.872**
	(0.429)
resid4	0.452
	(0.403)
propertyA122	0.446
	(0.325)
propertyA123	0.139
	(0.312)
propertyA124	0.950*
	(0.554)
age	-0.019
	(0.012)
other_installA142	-0.682
	(0.651)
other_installA143	-0.520
	(0.318)
housingA152	-0.391
	(0.314)
housingA153	-0.221
	(0.636)
other_credits2	0.455
	(0.312)
other_credits3	0.073
	(0.848)
other_credits4	0.851
	(1.384)
jobA172	-0.169
	(0.807)
jobA173	-0.198
	(0.779)
jobA174	0.257
	(0.807)
num_depend2	0.261
	(0.328)
telephA192	-0.372
	(0.265)
foreignA202	-2.272***
	(0.854)
Constant	-0.051
	(1.309)
Observations	600
Log Likelihood	-270.011
Akaike Inf. Crit.	650.022
Note:	*p<0.1; **p<0.05; ***p<0.01

>

Nyt voimme arvioida mallin ennustustarkkuutta laskemalla ROC-käyrän.
R-koodi:

testdata$model1p <- predict(model1, testdata, type = "response") # predict-komento laskee mallilla ennusteet testdatan arvoilla ja muodostaa testdataan uuden "model1p"-sarakkeen
library(ROCR)
model1points <- prediction(testdata$model1p, testdata$default) #vertaa ennusteita oikeisiin arvoihin
performancemodel1 <- performance(model1points,"tpr","fpr")
plot(performancemodel1, col = "blue")     #kuvaa ROC-käyrän
abline(0,1,lty=5)                                         #tekee 45 asteen suoran viivan

ROC-käyrän muotoja ja mallin onnistumistarkkuutta havainnollistaa alla oleva kuva.

Lähde: http://gim.unmc.edu/dxtests/roccomp.jpg

ROC-käyrässä True positive rate (TPR) tarkoittaa lukua, joka saadaan, kun jaetaan onnistuneet luotonhylkäykset, kun luotonsaaja on maksukyvytön, kaikilla maksukyvyttömien määrällä. False positive rate (FPR) tarkoittaa lukua, joka saadaan, kun jaetaan luoton hylkäykset, kun asiakas olisi ollut maksukykyinen, kaikkien maksukykyisten määrällä. Käyrän arvot saadaan, kun käytetään erilaisia ns. cut-off -arvoja ja lasketaan TPR ja FPR näillä arvoilla. Esim. kun cut-off-arvo on 0.1, niin luotot hylätään, jos asiakkaan maksukyvyttömyyden todennäköisyys on yli tämän. Hylätyistä luotoista lasketaan sitten oikein ja väärin luokitellut, jolloin saadaan ROC-käyrän piste tälle cut-off-arvolle.  ROC-käyrästä näemme, että cut-off arvo alkaa luvusta 1 ja oikealle mentäessä se laskee aina nollaan saakka, jolloin käyrät kohtaavat. Lasken vielä AUC-arvon, joka kertoo, kuinka paljon aluetta ROC-käyrän alle jää. Tämä auttaa ROC-käyrien vertailussa, kun ero ei visuaalisesti ole aivan selvä. ROC-käyrille saa laskettua myös luottamusvälit. 
R-koodi:

model1_auc <- performance(model1points, "auc")

Mallin AUC-arvo on 0.7898

Lasken myös Kolmogorov-Smirnov (KS) -arvon
R-koodi:


ks1 <- max(attr(performancemodel1, "y.values")[[1]] - (attr(performancemodel1 "x.values")[[1]]))



KS-arvo kertoo FPR:n ja TPR:n kumulatiivisen erotuksen maksimikohdan. KS-arvo 1 tarkoittaa, että malli on onnistunut jakamaan maksukyvyttömät ja maksukykyiset selkeästi kahteen eri luokkaan. KS-arvo 0 tarkoittaisi, ettei malli ole onnistunut luokkiin jakamisessa. KS-arvo 1 voitaisiin saada esim. jos kaikki maksukykyiset saisivat arvon <0.5 ja kaikki maksukyvyttömät arvon >= 0.5. Tällainen KS-arvo tarkoittaisi, että voisimme laittaa cut-off arvoksi 0.5, jolloin hylkäisimme kaikki maksukyvyttömien lainat ja myöntäisimme kaikki maksukykyisten lainat. 

KS-arvosta löytyy lisää tietoa mm. täältä. 

Kokeilin myös mallia, jossa jatkuvat muuttujat (age, duration ja amount) oli muutettu kategorisiksi, mutta mallin tulokset eivät eronneet merkittävästi. Mallia sai paranneltua, kun lisäsi siihen interaction-muuttujia. Potentiaalisia interaction-muuttujia löytää Machine learningiin (Koneoppiminen) kuuluvan Random Forestin Variable importance -listasta. Parhaan mallin löysin käyttämällä Koneoppimisen menetelmiä, joilla sain AUC-arvon, joka oli suurempi kuin D. Sharman paperissaan saama suurin arvo. Koneoppiminen on erittäin hyödyllinen menetelmä mallintamisessa, koska sen keinoin saadaan useista malleista luotua ennusteita, jotka peittoavat yksittäisten mallien tulokset. Näistä menetelmistä kirjoittelen jatkossa tarkemmin. 

Aikasarjan autokorrelaatio

Autokorrelaatiokuvat ACF ja PACF ovat olennaisessa osassa, kun tutkitaan aikasarjoja ja halutaan ennustaa. On siis tärkeää ymmärtää, miten kuvissa esiintyvät luvut syntyvät. Katsotaan nyt myyntien aikasarjaa, jota käytin SARIMA-esimerkissä ja tutkitaan autokorrelaatiota tarkemmin. Olen lisännyt kuvaan viivan, joka kuvaa myyntien keskiarvoa.

Alla myös aikasarjan autokorrelaatiokuvat.

Autokorrelaatio lasketaan alla olevalla kaavalla. 

Autokorrelaation suuruuteen vaikuttaa etenkin y(t):n ja y(t-1) välinen kovarianssi. Ensimmäistä kuvaa katsomalla näemme, että noin puolet arvoista on keskiarvon alapuolella ja noin puolet keskiarvon yläpuolella. Tästä johtuen, kun laskemme y(t):n autokorrelaatiota ensimmäiselle viiveelle, niin y(t):n ja y(t-1):n kovarianssi on jatkuvasti positiivinen. y(t):n varianssi on luonnollisesti positiivinen, joten tällöin myös autokorrelaatio on ensimmäiselle viiveelle positiivinen. 

Annan esimerkin autokorrelaation laskemisesta ensimmäiselle viiveelle. Ajatellaan, että meillä on aikasarja, joka muodostuu arvoista 1,2,3 ja 4, ja aikaväli on 01/2015-04/2015 (eli tammikuun arvo on 1, helmikuun arvo on 2, maaliskuun arvo on 3 ja huhtikuun arvo on 4). Olemme kiinnostuneita siitä korreloiko aikasarja ensimmäisen viiveensä kanssa, eli selittääkö tammikuun arvo helmikuun arvoa ja helmikuun arvo maaliskuun arvoa ja maaliskuun arvo huhtikuun arvoa. Lasketaan ensin kovarianssi y(t:):lle ja y(t-1):lle. Aikasarjan keskiarvo on (1+2+3+4)/4= 2.5 joten voimme kirjottaa kovarianssin seuraavasti: (2-2.5)(1-2.5)+(3-2.5)(2-2.5)+(4-2.5)(3-2.5)= 1.25

Lasketaan myös varianssi: (1-2.5)^2 + (2-2.5)^2 + (3-2.5)^2 + (4-2.5)^2 = 5

Autokorrelaatio ensimmäiselle viiveelle on siis 1.25/5=0.25

Lasketaan autokorrelaatiot myös viiveille 2 ja 3. Ensimmäiseksi lasketaan kovarianssi 2. viiveelle (huom.varianssi pysyy samana viiveistä riippumatta): (3-2.5)(1-2.5)+(4-2.5)(2-2.5)= -1.5 eli autokorrelaatio 2. viiveelle on: -1.5/5= -0.3

Ja lopuksi 3. viiveen kovarianssi: (4-2.5)(1-2.5)= -2.25 eli autokorrelaatio 3. viiveelle on: -2.25/5= -0.45

Positiivinen autokorrelaatio tarkoittaa sitä, että havainnot valitulla viivellä ovat suureksi osaksi saman merkkisiä, eli joko keskiarvoa pienempiä tai suurempia. Autokorrelaatiota aiheuttavat erityisesti mallista puuttuvat muuttujat sekä vääränlaisen mallin käyttö. Puuttuva muuttuja voi olla esim. AR(1) termi, jota emme ole vielä lisänneet malliin. Vääranlaisella mallillla tarkoitan tässä mallia, jossa on vielä yksikköjuuri. Kuten SARIMA-esimerkissä näimme, niin aikasarjassa on ennen kausidifferensointia vahva autokorrelaatio, joka vähenee huomattavasti tämän toimenpiteen jälkeen. Toisaalta kaikki autokorrelaatio ei katoa mallista ja tarvitsemme puuttuvat muuttujat malliin. On tärkeää tiedostaa, että ARIMA-ennustaminen toimii, koska aikasarjoissa eri havainnot eivät usein ole täysin itsenäisiä. Eli aiemman kuukauden havainto tarjoaa usein informaatiota tulevan kuukauden havainnosta. Itseasiassa "weakly" stationaarisessa AR(1) mallissa ensimmäisen viiveen autokorrelaatio on sama kuin mallissa esiintyvän y(t-1) kerroin, eli ensimmäinen viive vaikuttaa y(t):n aina samalla tavalla. Jos havainnot olisivat täysin itsenäisiä ja riippumattomia toisistaan, eli ns.white noise, niin emme voisi ennustaa kuin, että aikasarja pysyy keskiarvossaan. (White noisella tarkoitetaan sarjaa, jonka havainnot pyörivät keskiarvonsa ympärillä täysin sattumanvaraisesti, joten aikaisemmat arvot eivät tarjoa informaatiota tulevasta). Tämä on hyvä pitää mielessä, jos haluaa ennustaa esim. osakemarkkinoita. Tarjoavatko osakkeen historialliset arvostukset oikeasti relevanttia tietoa tulevista kursseista vai onko paras ennustus vain viimeisin kurssi?

PACF

Osittaisautokorrelaatio mittaa lineaarista riippuvuutta tietyllä viiveellä, kun aikaisempien viiveiden vaikutus on otettu pois. Tästä johtuen aikasarjan autokorrelaatio ensimmäisellä viiveellä on sama kuin osittaisautokorrelaatio ensimmäisellä viiveellä. Osittaisautokorrelaatiokuvat kertovat kannattaako AR-termejä laittaa ARIMA-/SARIMA -malleihin. Osittaisautokorrelaatio siis kertoo, onko esim. aikasarjan 3.viiveellä (y(t-3) lineaarista riippuvuutta y(t):n kanssa. Jos näin ei ole, niin aikasarjamalliin ei kannata laittaa AR(3) termiä.

Osittaisautokorrelaatio kuvataan usein, niin että osittaisiautokorrelaatio viiveellä h olisi sama asia kuin kyseisen viiveen regressiokerroin mallissa, jossa ovat sen lisäksi sitä edeltävät viiveet. AR-malli olisi siis viiveen h=3 tapauksessa seuraava: y(t)=b0+b1y(t-1)+b2y(t-2)+b3y(t-3) ja osittaisautokorrelaatio viiveelle 3 olisi kerroin b3.  Vastaavasti 4.viiveen AR-malli olisi y(t)=b0+b1y(t-1)+b2y(t-2)+b3y(t-3)+b4y(t-4) ja sen osittaisautokorrelaatio olisi kerroin b4. Tätä menetelmää ei kuitenkaan R:ssä käytetä, vaan kyseessä on Durbin-Levinson Rekursio, joka on algoritmi, jolla vastaavat kertoimet (osittaisautokorrelaatiot) lasketaan. Tästä laskentatavasta voit lukea lisää esim. täältä

Seuraavaksi simuloin R:llä ei-stationaarisen aikasarjan, jotta näemme, kuinka perinteisen AR-mallin kertoimet ja R:n osittaisautokorrelaatioluvut eroavat toisistaan huomattavasti. Simuloin myöhemmin myös stationaarisen aikasarjan ja näemme, että AR-mallin kertoimet ovat hyvin lähellä osittaisiautokorrelaatiolukuja.

R-koodi ei- stationaariselle aikasarjalle:

set.seed(1)                                                    #set.seed(1):lla saat samat tulokset kuin tässä
x <- w <- rnorm(500)
for (t in 2:500) x[t]<-x[t-1]+w[t]                 # Tehdään random walk-aikasarja
par(mfrow=c(2,1))                                       #saadaan kaksi kuvaa samaan ruutuun
plot(x,type="l")
Pacf(x,lag.max=3)

Pacf(x,plot=FALSE,lag.max=3)                #saadaan tulostettua arvot kolmelle viiveelle.

Kuva (ylempi) ei-stationaarinen aikasarja. Kuva (alempi) Osittaisautokorrelaatio

Partial autocorrelations of series ‘x’, by lag

     1        2          3
 0.972     0.009   -0.023

Yllä osittaisautokorrelaatioluvut. Nyt lasken R:n Arima-funktiota käyttäen vastaaville viiveille arviot, eli teen AR-mallit 1.viiveelle, 2.viiveelle ja 3.viiveelle.
R-koodi:

Arima(x,order=c(1,0,0))                                       #Saadaan 1.viiveen osittaisautokorrelaatio
Arima(x,order=c(2,0,0))                                       #Saadaan 2.viiveen osittaisautokorrelaatio
Arima(x,order=c(3,0,0))                                       #Saadaan 3. viiveen osittaisautokorrelaatio

Näin saadut osittaisautokorrelaatiot ovat

1          2           3
0.979   0.013   -0.007

Näemme, etteivät arviot ole kovin lähellä toisiaan.

Nyt muodostan stationaarisen aikasarjan ja teen samat toiminnot kuin ei-stationaariselle, eli kuva aikasarjasta ja osittaisautokorrelaatiosta sekä osittaisautokorrelaatioluvut.
R-koodi:

set.seed(1)
w<-rnorm(500)
par(mfrow=c(2,1))
plot(w,type="l")
Pacf(w,lag.max=3)
Pacf(w,plot=FALSE,lag.max=3)

Kuva (ylempi) stationaarinen aikasarja. Kuva (alempi) Osittaisautokorrelaatio

Saadut osittaisautokorrelaatioluvut ovat:

     1         2           3
-0.027    -0.004   -0.040                                                              

Nyt muodostan Arima-funktiota käyttäen AR-mallit, kuten aiemmin. 

R-koodi:

Arima(x,order=c(1,0,0))

Arima(x,order=c(2,0,0))

Arima(x,order=c(3,0,0))

Nyt saadut osittaisautokorrelaatioluvut ovat:

1            2             3

-0.027  -0.004     -0,040

Luvut ovat siis pyöristettyinä samat. Stationaarisessa aikasarjassa osittaisautokorrelaatiot ovat siis lähelle AR-mallin ennustamia, mutta eivät aivan samoja. Esim. 5.viiveelle lasketut arvot ovat 0.022 (Pacf) ja 0.026 (AR(5)). Intuitiivisesti osittaisautokorrelaation voi kuitenkin ymmärtää juuri AR-mallin tavoin. Ajatellaan esimerkiksi kahden mallin eroa, kun haluamme laskea 2.viiveen autokorrelaation. 1. malli: y(t)=y(t-2) ja 2. malli y(t)=y(t-1)+y(t-2)

1.mallissa emme ota huomioon 1.viiveen vaikutusta, kun taas toisessa otamme. Koska haluamme tietää vain 2.viiveen vaikutuksen y(t):hen, pitää y(t-1) ottaa mukaan, jotta sen vaikutus saadaan myös huomioitua.
























Kirjoittaja on rekrytointisivusto: www.dataverkosto.com perustaja. 

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

Tässä blogissa käsittelen kvantitatiivisia tutkimusmenetelmiä, kuten ARIMA:a, SARIMA:a ja Dynaamista regressiomallia. Käytän mm. SARIMA:a mainoskampanjan vaikutuksen analysointiin, jolloin tavoitteenani on pyrkiä vastaamaan kysymykseen "Onko mainoskampanjalla saatu aikaan lisää myyntiä?". Aineistona käytän tällöin aikasarjoja, eli aineistoja, jotka ovat ajallisessa järjetyksessä. Aikasarja-aineiston käyttö mahdollistaa kampanjan tutkimisen tietyllä aikavälillä. SARIMA mahdollistaa kampanjan vaikutuksen, kun saatavilla ovat ainoastaan tuotteen historialliset myyntimäärät, eli esim. Tuotteen X myynti aikavälillä 2005-2014.  Esittelen myöhemmin myös Dynaamista regressiomallia hyödyntävän tavan, jolla voimme tutkia kampanjan vaikutusta, kun myös muuta informaatiota (esim. mainostuskulut) on saatavilla. Pyrin esittämään mallit ymmärrettävästi ja yksinkertaisesti, koska olen huomannut, ettei tällaista tekstiä suomenkielellä vielä ole. Tämän vuoksi en pureudu mallien toimintaan juurikaan teoriatasolla, mutta selitän mallien toiminnat sanallisesti, jotta lukija saa intuitiivisen kuvan niiden toimintaperiaatteista. Aineiston käsittelyyn jne. käytän R-ohjelmaa.

Ensimmäisenä esittelen ARIMA-mallin. ARIMA:a käytetään aikasarja-analyysiin ja se soveltuu ennustamiseen. ARIMA perustuu siihen, että tietyn aikasarjan kehitystä (esim. myynnit 2005-2013) selitetään sen omalla historialla, eli voimme ennustaa tuotteen ensi kuukauden myyntiä tämän kuukauden myynneillä. ARIMA koostuu AR (Autoregressive)- ja MA (Moving average)  -prosesseista, jotka vaikuttavat käytännössä siihen, kuinka malli reagoi omiin historiallisiin arvoihinsa. Puhdas AR-prosessi tarkoittaa, että myynnit tänä kuukautena ovat riippuvaisia ainoastaan aiempien kuukausien myynneistä. AR-prosessit merkitään muodossa AR(x), jossa x tarkoittaa sitä, kuinka monta aiempaa myyntihistoriaa käytetään tulevan arvon ennustamiseen. Esim. AR(1) tarkoittaa, että tulevan kuukauden myyntiä selitetään vain tämän kuukauden myynnillä. AR(2) tarkoittaisi taas, että myös viime kuun myynti vaikuttaisi ennustettuun arvoon.

MA-prosessi on prosessi, joka on mielestäni hieman hankalampi ymmärtää intuitiivisesti. ARIMA-mallissa se tarkoittaa yleisesti shokkia, joka kohdistuu esim. myyntiin, mutta joka ei selity AR(x) -prosessilla. Esimerkiksi voimme ajatella kenkien myyntiä, jossa tapahtuu positiivinen kysyntäshokki esim. mainoskampanjan takia. Mikäli meillä olisi puhdas AR(1) -malli ja selittäisimme kenkien myyntiä ainoastaan viimeisimmän kuukauden myynneillä, niin kysyntäshokki jäisi vähälle huomiolle. AR(1) malli ei ota shokkia erikseen huomioon, kun se ennustaa seuraavaa kuukautta. Se on vain riippuvainen kenkien myynnistä, eli koska shokki lisäsi kenkien myyntiä, niin AR(1) ennustaa muutosta seuraavalla kuukaudelle. Lisäämällä MA -termin malliimme voimme hienosäätää mallia ja näin ottaa shokin erikseen huomioon, kun teemme ennustetta. Lisäämällä MA(1) -termin AR(1)-malliimme voimme ennustaa shokin vaikuttavan ensi kuukauden myyntiin. Lisäämällä MA(2) -termin malliimme voimme ennustaa shokin vaikuttavan vielä kahden kuukauden päähän. MA-prosessilla saamme siis normaalin AR-mallin sijaan mahdollisesti tarkempia ennustuksia. MA-prosessin avulla voimme esim. ennustaa, että myynnit laskevat shokin jälkeen nopeammin kuin pelkkä AR(1) -malli antaisi ymmärtää.

On huomattava, että MA-termin vaikutus on yleensä lyhytaikainen ennustuksissa, koska sitä varten tarvitsemme mallin ennustaman arvon sekä oikean arvon. Kun käytämme ARIMA:, tarvitsemme oikeaa historiallista dataa, josta ensin muodostamme mallimme. ARIMA-malli muodostetaan automaattisesti R-ohjelman toiminnoilla mallintamaan oikea dataa. Käytännössä ohjelma laskee AR:n ja MA:n parametrien arvot, jotka minimoivat mallin ja oikeiden arvojen erotuksen. Jos haluamme ennustaa ARIMA-mallilla esim.3 kuukautta eteenpäin ja meillä on AR(1)- ja MA(1) -termit mallissa, MA(1)-termi vaikuttaa ainoastaan ensimmäiseen ennusteeseen, eli aina kun meillä on MA mukana, niin tarvitsemme oikean historallisen arvon. AR(1)-prosessi taas pystyy jatkumaan, vaikka meillä ei ole historiallista arvoa, koska se saa arvonsa aiemmasta mallin arvosta. Esim. jos AR(1) -malli ennustaa helmikuulle 2009 myynniksi 150000€, niin sen ennuste maaliskuulle perustuu puhtaasti tähän jo ennustettuun arvoon, joten pitempien ennusteiden kanssa pitää olla varovainen.

ARIMA-esimerkin jälkeen käsittelen SARIMA:a, joka eroaa ARIMA:sta, koska kausivaihtelu (Seasonality) on otettu huomioon. Alla olevassa kuvassa näkyy kausivaihtelu, eli vaihtelu, joka on joka vuosi lähes samanlainen. Esim. kaupan alalla joulukuiset myyntipiikit ovat yleisiä.

Kuva 1. Firman X Oy myynnit 1999-2010.

Kuvasta näkyy myös vuosittaisen kausivaihtelun lisäksi selkeä nouseva trendi, joka tulee sekin ottaa huomioon ARIMA:ssa ja SARIMA:ssa. Kausivaihtelu ja trendi tulee ottaa huomioon, jotta saamme mallimme käyttäytymään, kuten historiallisesti myynnit kehittyvät. Jos emme ottaisi näitä huomioon, niin ennusteemme yllä olevaan kuvaan ei muistuttaisi historiallista muotoa, vaan olisi nouseva/laskeva/tasainen viiva. Toisaalta trendiä ei ole pakko ottaa ennusteissa huomioon, jos uskoo, ettei esim. kuvan nouseva trendi jatku. Mikäli trendi lisätään malliin, niin kuvan tapauksessa myös ennusteemme tulee nousemaan. Trendin lisäys voisi siis johtaa liian optimistisiin myyntiodotuksiin. ARIMA/SARIMA ei ota kantaa siihen, miksi kausivaihtelut ovat mitä ovat tai miksi trendi kehittyy tietyllä tavalla. Voimme kuitenkin hyödyntää kuvan informaatiota ja olettaa, että kausivaihtelut noudattavat samaa muotoa myös tulevaisuudessa. Jos emme tekisi näin ja ennustaisimme joulukuun myyntiä, niin arviomme menisi luultavasti alakanttiin. Ottamalla kausivaihtelun huomioon saamme ARIMA/SARIMA-mallimme korjaamaan ennusteen joulukuulle.

Käsittelen myös Dynaamista regressiomallia. Erityisesti esittelen tapauksen, jossa regressiomallin virhetermi on ARIMA-prosessi. Tällöin tavallinen regressiomalli ei riitä, koska virhetermeille asetetut oletukset eivät päde, joten mallista saadut tulokset eivät ole luotettavia. Dynaaminen regressiomalli tarkoittaa yleisesti sitä, että regressiomalli sisältää selittävänä tekijänä selitettävän tekijän viivästetyn arvon (lag:n), eli esim. myyntejä selittää myös aiemman kuukauden arvo. Aivan kuten ARIMA/SARIMA-malleissa. Dynaaminen regressiomalli antaa meille mahdollisuuden lisätä ARIMA/SARIMA-prosessiin myös muita selittäviä tekijöitä, eli voimme lisätä malliin esim. kilpailijan hinnoittelumuutokset ja omat markkinointikampanjat.