R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Aikasarjaennusteiden laatiminen

Korkeasaari on eläintarha, joka sijaitsee Helsingissä. Korkeasaareen pääsee joko vesibussilla tai maitse Mustikkamaan kautta. Eläintarha on suosittu vierailukohde, erityisesti kesäaikaan. Korkeasaaren kävijämäärät on julkaistu avoimena datana CC-BY-4 -lisenssillä.

Aineistossa on pilkottu kävijämäärät vuosille 2010-2014 kuukausittain ja saapumisreiteittäin sekä lipputyypeittäin. Tästä muodostuu varsin monipuolinen aikasarja-aineisto, joka soveltuu myös erilaisten aikasarjaennusteiden tekemiseen. Rob J Hyndmanin ja George Athanasopoulosin kirja Forecasting: Principles and Practise on tässä suureksi avuksi.

Työssä tarvitaan seuraavat R-paketit:

# Varsinaiset analyysit
library(xlsx)
library(fpp2)
library(lubridate)
library(ForecastComb)

Aineisto

Eri vuosien ja kuukausien päiväkohtaiset kävijämäärät pötkötettiin yhteen samaan tiedostoon, yhteenveto.xlsx. Tämän jälkeen laskettiin aineistosta viikkokohtaiset kävijämäärät:

# Luetaan data
d <- read.table("C:\\Users\\lenovo\\Desktop\\korkeasaari\\yhteenveto.xlsx", sheetName="Taul2")

# Muokataan dataa
d$date <- as.Date(d$pvm, format="%d.%m.%Y")
d$week <- isoweek(d$date)
d$year <- year(d$date)
d <- na.omit(d)
d$week[d$week==53] <- 52
a <- aggregate(d$kaikki_yht, list(d$week, d$year), sum)

Viikkotason kävijämääristä muodostetaan aikasarja seuraavasti. Vuoden 2016 havainnot jätetään erikseen testiaineistoksi, jolla voidaan arvioida ennusteiden osuvuutta.

# Tehdään kokonaismäärästä aikasarja
# Jätetään vuosi 2016 ennustetarkkuuden mittaamiselle erikseen
y <- ts(a$x[1:c(52*6)], start=c(2010, 1), frequency=52)

Aikasarjasta voidaan tehdä erilaisia kaavioita, jotka esittävät koko aikasarjan (autoplot) tai niin sanottu vuodenaikaisplotti (ggseasonplot), jossa eri vuosien havainnot on esitetty samassa kuvassa:

 
# Perusplotteja
autoplot(y)
ggseasonplot(y)
ggsubseriesplot(y)
ggAcf(y)

Tuloksena syntyvät seuraavat kuvat:

Kuvista on helppo havaita, että kävijämäärät ovat suurimmillaan kesällä.

Ennusteet

Laaditaan aikasarjaennusteet neljällä eri menetelmällä, ja plotataan niiden ennusteet vuodelle 2015 samaan kuvaan oikeiden vuoden 2015 havaintojen kanssa. Oikeat havainnot on esitetty kuvassa mustalla viivalle.

# Naiivi
fit0 <- snaive(y, 52)
plot(x=1:52, y=a$x[261:nrow(a)], type="l", xlab="Viikko", ylab="Kävijämäärä")
lines(x=1:52, y=as.data.frame(snaive(y, 52))[,1], type="l", col="red")

# auto.arima
arimaMod <- auto.arima(y, stepwise=FALSE, approximation=FALSE)
arimaMod.Fr <-forecast(arimaMod,h=52) 
lines(x=1:52, y=arimaMod.Fr$mean, col="blue")

# decomposition
fcast <- stlf(y, method='naive')
lines(x=1:52, y=fcast$mean[1:52], col="green")

# neural network
fitn <- nnetar(obs)
lines(x=1:52, y=forecast(fitn, 52)$mean, col="grey50")

Tuloksena on seuraava kuva:

Yhden ennusteen (tässä arima) voi myös plotata luottamusväleineen esimerkiksi seuraavalla komennolla

plot(arimaMod.Fr)

joka muodostaa seuraavan kuvan:

Ennusteiden yhdistäminen

Yksittäiset aikasarjaennusteet voidaan myös yhdistää, vähän samaan tapaan kuin ensemble-malleissa tehdään. Paketti forecastComb tarjoaa tähän työkaluja.

Tehdään yhdistely kahdella eri menetelmällä:

obs <- a$x[1:260]
preds <- matrix(c(fitted(fit0), fitted(arimaMod.Fr), fitted(fcast), fitted(fitn)), ncol=4, nrow=260)
test_obs <- a$x[261:312]
test_preds <- matrix(c(fit0$mean, arimaMod.Fr$mean, fcast$mean[1:52], forecast(fitn, 52)$mean), ncol=4, nrow=52)
fc <- foreccomb(obs, preds)
fc$Forecasts_Train <- as.ts(fc$Forecasts_Train)
resNG <- comb_NG(fc)
resBG <- comb_BG(fc)
respNG <- predict(resNG, test_preds)
respBG <- predict(resBG, test_preds)

Pyöräytetään lopuksi vielä kaikki ennusteet samaan taulukkoon ja piirretään siitä viivakaavio:

vertailu <- data.frame(obs=a$x[261:312], snaive=fit0$mean, arima=arimaMod.Fr$mean, stlf=fcast$mean[1:52], nn=forecast(fitn, 52)$mean, NG=respNG, BG=respBG)
matplot(vertailu, type="l", lty=1)

Tuloksena oleva kuva osoittaa jälleen eri ennusteiden olevan varsin lähellä toisiaan:

Mikä ennusteista sitten on kaikkein lähimpänä havaittuja arvoja testidatassa (vuosi 2015)?

Ennusteen osuvuus

Ennusteiden osuvuutta voidaan tutkia esimerkiksi komennolla accuracy(). Se voidaan helposti ajaa kaikille eri malleille yllä muodostetusta taulukosta seuraavasti:

acc <- apply(vertailu[,-1], 2, function(x){ accuracy(vertailu[,1], x)} )
rownames(acc) <- c("ME", "RMSE", "MAE", "MPE", "MAPE")

Tuloksen muodostuvassa taulukossa kukin rivi on yksi tarkkuuden mittaamiseen tarkoitettu suure.

         snaive      arima       stlf         nn          NG          BG
ME   -734.65385 -250.74776 -1405.2343 -317.86221 -501.081184 -709.199801
RMSE 3639.84851 3144.43747  3511.3019 3154.40535 3098.004287 3179.234575
MAE  2401.96154 2147.39850  2266.3893 2376.21863 2050.404338 2029.617821
MPE   -11.60592    8.06956  -384.8611   11.48996    3.617586   -3.745891
MAPE   38.47657   33.51295   512.2381   37.52990   32.651673   30.783737

Riippuen tarkkuuden mittaamiseen käytetystä suureesta, eri mallit toimivat parhaiten. Esimerkiksi arima ja neuroverkko vaikuttavat toimivan parhaiten, jos virheen suuruutta mitataan keskimääräisellä virheen suuruudella (ME eli tässä sum(pred - obs) / 52). Jos taas virheen etumerkkiä ei huomioida, jolloin positiiviset ja negatiiviset virheet eivät kumoa toisiaan (MAE eli tässä sum(abs(pred - obs)) / 52)), yhdistelmämenetelmät NG ja BG tuntuvat toimivan paremmin. Riippuu sitten käyttötarkoituksesta, mikä tarkkuuden arviointiin käytetty menetelmä soveltuu parhaiten.

Hierarkkiset aikasarjat

Yllä on ennustettu vain kävijöiden kokonaismäärää. Aineistossa on kuitenkin saatavilla myös eri saapumisreittien ja lipputyyppien mukaan jaoteltu kävijämäärä. Esimerkiksi maitse saapuneet pääsynsä maksaneet, perhelippulaiset ja ilmaislippulaiset summautuvat yhteen maitse saapuneisiin, joka sitten tietysti edelleen summautuu osana kävijöiden kokonaismäärään. Tämä hierarkkinen rakenne voidaan ottaa huomioon ennusteita tehtäessä.

Hierarkkisten aikasarja-aineiston käsittelyyn löytyy R-paketti hts, ja lisäksi aiemmin mainitun kirjan 10. lukuh käsittelee nimenomaan tätä asiaa.

Saattaa olla, että innostun vielä kokeilemaan hierarkkisen aikasarjan perusteella ennusteiden tekemistä tällä samalla aineistolla, mutta jollen, niin lisätietoa löytyy siis yllä mainituista lähteistä.

Yhteenveto

Aikasarjaennusteita voidaan tehdä useilla erilaisilla menetelmillä. Eri menetelmien ennusteiden yhdistäminen saattaa parantaa kokonaisennustetta pienentämällä virheen suuruutta. R-paketit fpp2 ja forecastComb tarjoavat helppokäyttöisiä työkaluja näihin tilanteisiin.