## Tarvitset paketin  forecast (Packages->install package(s) -> "paikka" --> forecast) 

## Haetaan aineisto  - tässä tapauksessa Lake (voit olla tallentanut haluamasi tiedoston minne vain)
## header=true olettaa, että tiedoston ensimmäisellä rivillä on muuttujien nimet

data=read.table(file.choose(), header=TRUE)  ## Tämän komennon pitäisi antaa sinun valita kansioistasi haluamasi tiedosto

## Vaihtoehtoinen tapa. Muuta tiedoston polku oikeaksi ja poista kommentointi jos haluat käyttää tätä menetelmää.
##
## data=read.table("C:/Opiskelu/Tilastotiede/Stationaariset 2015/Laskuharjoitukset/Lake.txt",header=TRUE)

y = ts(data$Y)  ## Muodostetaan aikasarja. Komento etsii sarakkeen Y aineistosta datan, muuttaa sen aikasarjaksi ja nimeää y:ksi.

## Vaihtoehtoinen tapa, kun halutaan antaa ensimmäisen havainnon alkamisajankohta (start=1959) 
## ja havaintojen lkm/aikayksikkö (freq=4). Ks. ?start

y=ts(data$Y, start=1892,freq=1) ## Aineiston alkuajankohta ja frekvenssi on ilmoitettu. freq = 1 tarkoittaa vuosihavaintoja,
			    ## freq = 4 neljännesvuosihavaintoja, freq = 12 kuukausihavaintoja, jne ...
			    ## Jos aikayksikköä kohti on enemmän kuin yksi havainto, täytyy alkuajankohta ilmoittaa vektorina. 
			   ## Syyskuu 1959 on start = c(1959, 3).


## Aikasarjan piirtäminen (ks. R-koodi_1)
## Muunnokset tarvittaessa (ks. R-koodi_1) 


## MALLIN VALINTA SU-ESTIMOINTIIN PERUSTUVIEN MALLINVALINTAKRITEERIEN AVULLA

## EI-KAUSIVAIHTELUTAPAUS

## ARMA(p,q)-mallin asteiden valinnassa voidaan käyttää auto.arima-koodia, joka perustuu 
## mallinvalintakriteerien käyttöön suurimman uskottavuuden estimointimenetelmässä. Ks. ?auto.arima 

## Jos diffrensointi tarpeen, se oletetaan seuraavssa tehtävän 'manuaalisesti' etukäteen, jolloin asteet p ja
## q valitaan differensoidulle sarjalle (eli auto.arima-koodin mahdollisuutta valita differensointiastekin
## ei käytetä). Ks. ?diff

## Edellä todetun mukaisesti, kiinnitetään alempana d=0
## p = AR-viipymien määrä, q=MA-viipymien määrä
## 'max.p'  ja 'max.q'  valitsevat AR- ja MA-viipymien maksimimäärät; tässä molemmat 5
## 'max.order' on summan p+q  maksimiarvo; tässä 5
## 'seasonal = FALSE' rajoittaa mallinvalinnan ei-kausivaihtelu- eli ARMA(p,q)-malleihin (kausivaihtelumalleja
## tarkastellaan monisteen jaksossa 4.6 ja niissä tarvittava koodi esitetään erikseen alempana)

## 'trace=TRUE', jos luettelo tarkastelluista ARMA-malleista tulostetaan, muuten 'trace=FALSE'
## 'stationary=TRUE', koska rajoitutaan malleihin, jotka ovat stationaarisia mahdollisen differensoinnin jälkeen
## 'stepwise=FALSE', kun kaikki mahdolliset mallivaihtoehdot käydään läpi
## 'approximation=FALSE', jos käytetään tarkkaa SU-estimointia; 'approximation=TRUE', jos käytetään approksimatiivista  
## SU-estimointia, joka on nopeampi, kun havaintoja on paljon ja tarkasteltavien mallien määrä suuri.

fit0 <- auto.arima(y, d=0, max.p=5, max.q=5, max.order=5, seasonal=FALSE, trace=TRUE, 
              		stationary=TRUE, stepwise=FALSE, approximation=FALSE)
fit0

## Huomaa, että koodi käy läpi myös mallivaihtoehdot, joissa ei ole vakiotermiä. Jos tällainen malli tulee 
## valituksi, mutta haluat sisällyttää malliin vakion, valitse 'trace=TRUE', tarkastele tuloksia vakion 
## sisältävien mallien osalta ja valitse niiden perusteella sopivimpana pitämäsi malli.
## Vastaavasti, jos vakiollinen malli tulee valituksi, etkä halua sisällyttää malliin vakiota.

## KAUSIVAIHTELUTAPAUS

## ARMA(p,q)x(P,Q)_s-mallin asteet voidaan valita auto.arima-koodilla, joka perustuu mallinvalintakriteerien 
## käyttöön suurimman uskottavuuden estimointimenetelmässä. Ks. ?auto.arima 

## Jos diffrensointi tarpeen, se oletetaan tehtävän 'manuaalisesti' etukäteen, jolloin asteet p, q, P ja Q 
## valitaan differensoidulle sarjalle (eli auto.arima-koodin mahdollisuutta valita differensointiasteetkin 
## ei käytetä). Kausidifferensoinnin saa komennolla diff(y,lag=s), jossa s=kausivaihtelun pituus. Ks. ?diff

## Huomaa, että auto.arima- ja esimoinnissa käytettävä Arima-koodi olettavat, että s on sarjan frekvenssi 
## (freq; ks.  y=ts(data$Y, start=1892, freq=1) edellä)

## Edellä todetun mukaisesti, valitaan alempana d=0 ja D=0
## p ja q kuten aikaisemmin, P=AR-kausiviipymien määrä, q=MA-kausiviipymien määrä
## 'max.p', 'max.q', 'max.P' ja  'max.Q'  valitsevat AR-osien ja MA-osien viipymien 
## maksimimäärät; tässä p=q=5 ja P=Q=2
## 'max.order' on summan p+q+P+Q  maksimiarvo; tässä 5
## 'trace', 'stationary', 'stepwise' ja 'approximation'  kuten ei-kausivaihtelutapauksessa

##  HUOM: Samat huomautukset kuin ei-kausivaihtelutapauksessa pätevät tässäkin. 

fit0 <- auto.arima(y, d=0, D=0, max.p=5, max.q=5,  max.P=2, max.Q=2, max.order=5, 
     		   trace=TRUE, stationary=TRUE, stepwise=FALSE,  approximation=FALSE)
fit0


## PARAMETRIEN ESTIMOINTI SU-ESTIMOINTI

## EI-KAUSIVAIHTELUTAPAUS

## Arima-koodilla estimoidaan haluttu ARIMA(p,d,q)-malli; ks. ?Arima
## Toisin kuin mallinvalinnassa edellä, ei differensointia tehdä tässä 'manuaalisesti' etukäteen, 
## koska ennusteet halutaan laskea differensoimattomalle sarjalle (jos halutaan ennustaa 
## differenssejä, käytetään tässäkin differensoitua sarjaa).

## p ja q kuten edellä ja d differensointiaste
## 'include.mean=TRUE', kun tapauksessa d=0 mallissa on vakiotermi. 
## 'include.drift=TRUE', kun tapauksessa d=1 mallissa on vakiotermi, jolloin (ennusteessa) on deterministinen
## lineaarinen trendi (ks. HT 5.3 ja 5.4(iv)) . Tällöin 'include.mean=TRUE' on merkityksetön.
## Jos valitaan 'd=0', 'include.mean=TRUE'  ja 'include.drift=TRUE', estimoidaan ARMA(p,q)-malli, jossa on
## deterministinen lineaarinen trendi eli ns. trendistationaarinen malli (ks. HT 5.1 ja 5.4(v)).  Jos haluat valita mallin 
## asteita olettaen deterministisen lineaarisen trendin, kannattaa trendi poistaa ennen auto.arima-koodin käyttöä. 
## 'transform.pars=TRUE' pakottaa parametrit stationaarisuusalueelle

fit<-Arima(y, order = c(p,d,q), include.mean=TRUE, include.drift=TRUE, transform.pars=TRUE)  
fit

## KAUSIVAIHTELUTAPAUS

## Jos käytetään differensointia eli estimoidaan ARIMA(p,d,q)x(P,D,Q)_s-malli, menetellään seuraavasti.
## p, q, P, Q  ja d  kuten edellä, D=kausidifferensointiaste
## 'include.mean=TRUE', kun tapauksessa d=0 ja D=0 mallissa on vakiotermi. 
## 'include.drift=TRUE', jos differensoidussa mallissa (eli d=1 ja/tai D=1) on vakiotermi, jolloin (ennusteessa) on 
## deterministinen lineaarinen trendi. Tällöin 'include.mean=TRUE' on merkityksetön. 
## Jos valitaan 'd=0', 'D=0' ja 'include.drift=FALSE', estimoidaan ARMA(p,q)x(P,D,Q)_s-malli, jossa deterministinen 
## lineaarinen trendi eli ns. trendistationaarinen malli (ks. HT 5.4(iv)). Tällöin 'include.mean=TRUE' on merkityksetön.
## ## 'transform.pars=TRUE' pakottaa parametrit stationaarisuusalueelle

fit<-Arima(y,order = c(p,d,q), seasonal = c(P,D,Q), include.mean=TRUE, include.drift=FALSE,  transform.pars=TRUE) 
fit

## Tiedoksi, että kausivaihtelu ja ei-kausivaihtelutapaus voidaan myös korvata yhdellä koodilla. Ks. ?Arima 


## VALITUN MALLIN DIAGNOSTIIKKAA 

fit<-Arima(y, order = c(p,d,q), include.mean=TRUE, include.drift=FALSE, transform.pars=TRUE)  
fit
tsdiag(fit)  ## Estimoidun mallin diagnostiikkaa 

## tai vastaava kausivaihtelutapaus

fit<-Arima(y,order = c(p,d,q), seasonal = c(P,D,Q), include.mean=TRUE, include.drift=FALSE,  transform.pars=TRUE) 
fit
tsdiag(fit)  ## Estimoidun mallin diagnostiikkaa 

## Jatko sama sekä kausivaihtelu- että ei-kausivaihtelutapauksissa

## Vaihtoehtoisesti residuaalien tutkiminen "manuaalisesti" (suositeltavaa ainakin niiltä osin, jotka eivät sisälly edelliseen).
## Koska Ljung-Box -testissä käytetään nyt residuaalisarjaa, vähennetään khi^2-jakauman vapausasteluvusta
## ARMA-asteiden summa eli p+q, jonka seuraavassa on oletettu olevan 2 ("fitdf")

resid=fit$residuals  ## "residuals" on yksi komponentti objektissa "fit". Se sisältää mallin residuaalit.

acor=acf(ts(resid), lag=40, type = "correlation", main="Residuaalien empiirinen autokorrelatiofunktio")
Box.test(resid, lag = 10, type = c("Ljung-Box"), fitdf = 2)  ## ks. ?Box.test

## Vastaavasti neliöityjen residuaalien ja McLeod-Li -testin tapauksessa (nyt khi^2-jakauman 
## vapausasteluvusta ei vähennetä mallin ARMA-asteiden summaa).

acor=acf(ts(resid^2), lag=40, type = "correlation", main="Neliöityjen residuaalien empiirinen autokorrelatiofunktio")  
Box.test(resid^2, lag = 10, type = c("Ljung-Box"), fitdf = 0)

## Suoritetaan standardointi ja muodostetaan standardoitujen residuaalien aikasarjakuva, histogrammi ja qq-plot

sresid=resid/sqrt(fit$sigma2)
plot(sresid)
hist(sresid, main="Standardoitujen residuaalien histogrammi") 
qqnorm(sresid, main="Standardoitujen residuaalien qq-plot")  ## Normaalisuusoletuksen tutkimiseksi
qqline(sresid)


## ESTIMOIDULLA MALLIlla ENNUSTAMINEN

## Estimoidulla ARIMA-mallilla ennustetaan forecast-koodin avulla; ks. ?forecast.Arima
## Saatavat ennusteet poikkeavat tapauksessa q>0 monisteessa esitetyistä ennusteista, koska ne perustuvat 
## prosessin äärelliseenhistoriaan eli havaittuun aikasarjaan. Erot näiden eri tapojen välillä ovat kuitenkin 
## pienet, ellei havaintoja ole kovin vähän.
## Seuraava koodi toimii sekä kausivaihtelu- että ei-kausivaihtelutapauksissa.

## h = ennustehorisontin pituus
## 'level', antaa mahdollisuuden valita ennusteen luottamusrajat; tässä 68% ja 95%
## 'fan = TRUE' esittää luottamusrajat edellisestä poikkeavaa grafiikkaa käyttäen
 
fcast <- forecast(fit, h=10, level=c(68,95), fan=FALSE)
fcast # Tulostaa ennusteeet ja luottamusrajat
plot(fcast) # Ennustekuva valituin asetuksin