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

## Haetaan aineisto  - tässä tapauksessa US_GDP (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 2014/Laskuharjoitukset/US_GDP.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=c(1984,1)) 
## ja havaintojen lkm/aikayksikkö (freq=4). 
## Ks. ?start

y=ts(data$Y, start=c(1984,1),freq=4)   

## 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. Helmikuu 1980 olisi start = c(1980, 2).

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

## Parametrien estimointi ja mallinvalintakriteerit uskottavuusfunktion avulla ?arima 

## arima-koodilla estimoidaan haluttu ARMA(p,q)-malli 
## p=AR osan viipymien määrä, q=MA osan viipymien määrä

mean=mean(y)  ## Keskistys tarvittaessa
yk=y-mean

fit<-arima(yk,order = c(p,0,q),include.mean = FALSE)  ## 'include.mean = TRUE', jos keskistystä ei tehdä (jolloin yk -> y).
fit

## tai haluttu ARMA(p,q)x(P,Q)-kausivaihtelumalli 
## p ja q kuten edellä, P=AR-kausiviipymien määrä, Q=MA-kausiviipymien määrä
## s=kausivaihtelujakson pituus
## kausidifferensointi tarvittaessa; ?diff
## Katso ?arima.

fit<-arima(yk,order = c(p,0,q),seasonal = list(order = c(P,0,Q), period = s),include.mean = FALSE) 
fit

## Voit kokeilla eri asteita laskemalla informaatiokriteereiden arvoja (tulostus esim. Info_AIC) 
## Huomaa, että R:n edellä tulostamat informaatiokriteereiden arvot lasketaan eri tavalla kuin kurssilla ja alla.

n_parameters=length(fit$coef)  ## "length" antaa vektorin pituuden. arima -funktio estimoi mallin jolle on annettu nimi fit. 

## Tämä objekti sisältää komponentteja, joita "voi poimia dollarimerkillä". fit$coef poimii vektorin,
##  joka sisältää estimoidut parametrien arvot. Funktion ohjesivulta löytyy lista komponenteista.

T=length(y)

Info_AIC=log(fit$sigma2) + 2*(n_parameters)/T
Info_BIC=log(fit$sigma2) + (n_parameters)*log(T)/T
Info_HQ=log(fit$sigma2) + 2*(n_parameters)*log(log(T))/T


## Valitsemasi mallin diagnostiikkaa

fit<-arima(y,order = c(p,0,q),include.mean = TRUE)  ## tai vastaava kausivaihtelumalli
tsdiag(fit)  ## mallin diagnostiikkaa 

## Vaihtoehtoisesti residuaalien tutkinta "manuaalisesti".
## Koska Ljung-Box -testissä käytetään nyt residuaalisarjaa, vähennetään khi^2-jakauman vapausasteluvusta
## asteiden summa eli p+q (tai kausivaihtelumallissa p+q+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 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 histogrammi ja qq-plot

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

## Estimoidulla mallilla ennustaminen ?predict.Arima

## Huomaa, että ennusteet lasketaan harjoitustehtävän 5.2 huomautuksessa kuvatulla tavalla, joka 
## poikkeaa monisteessa esitetystä tavasta. Erot näiden eri tapojen välillä ovat kuitenkin pienet,
## ellei havaintoja ole kovin vähän.

## Tapaus I: Ennustaminen viimeisestä havainnosta eteenpäin

od <- options(digits = 4)  ## Valitaan (haluttaessa) ennusteiden keskivirheissä käytettyjen desimaalien määrä (tässä 4)
fcast <- predict(fit, n.ahead = 8, se.fit = TRUE) ## "n.ahead = 8" ilmoittaa laskettavien ennusteiden lukumäärän (tässä 8), "se.fit = TRUE", kun keskivirheet palautetaan
fcast

## Ennusteiden ja niiden luottamusvälien graafinen esittäminen

f <- fcast$pred  ## Ennusteet vektoriin f
f <- f + mean    ## Lisätään ennusteisiin keskiarvo (sivuutetaan, jos keskistystä ei ole tehty)
se <- fcast$se   ## Ennusteiden keskivirheet vektoriin se

lv_ala <- f - 1.96*se  ## Vektori ennusteiden (tässä 95%:n) luottamusvälin alarajalle
lv_yla <- f + 1.96*se  ## Vektori ennusteiden (tässä 95%:n) luottamusvälin ylärajalle

forecast <- ts(c(y, f[NA]), start=c(1984,1), freq=4)  ## Havaittu sarja ja ennuste ensin vektoriin forecast, jota jatketaan ennustehorisontin loppuun. 
plot(forecast, ylim=c(min(y, lv_ala, lv_yla),max(y, lv_ala, lv_yla)), lwd=1, lty=1,
    ylab="", main="US_GDP-sarja ja sen ennuste",xlab="")  ## Valitaan kuvan ylä- ja alaraja, jotta luottamusväli mahtuu kuvaan (voit joutua laajentamaan rajoja)

lines(f,lwd=2, col="blue", lty=5)  ## Piirretään ennuste sinisellä
lines(lv_ala, col="red", lwd=2, lty=4)  ## Piirretään ennusteen alaraja punaisella
lines(lv_yla, col="red", lwd=2, lty=4)  ## Piirretään ennusteen yläraja punaisella
abline(v=c(2013,1), lty=3, lwd=1)   ## Piirretään pystysuora viiva, joka erottaa havainnot ja ennusteet. Ajankohdan määrittävän vektorin c(2013,1) valinnasta, ks. ?start

## Edellä "lwd" määrittelee viivan paksuuden, "lty" viivan tyypin ja "col" värin. Seuraavat vaihtoehdot ovat (ainakin) käytettävissä
## lwd: line width, positiivinen kokonaisluku oletusarvona 1
## lty: line type, vaihtoehtoina 0=blank, 1=solid (oletusarvo), 2=dashed, 3=dotted, 4=dotdash, 5=longdash, 6=twodash
## col: vaihtoehtoina "black", "red", "green3",  "blue" , "cyan", "magenta", "yellow", "gray"   


## Tapaus II: Ennusteiden vertaaminen valittuun määrään havaitun sarjan viimeisiä havaintoja

y1 = ts(y,start=c(1984,1),end=c(2008,4), freq=4)  ## Valitun estimointiperiodin (tässä 1984I - 2008IV) havainnot vektoriin y1

mean=mean(y1) ## Keskistys tarvittaessa
yk=y1-mean

## Estimoidaan haluttu ARMA(p,q)-malli 

fit1<-arima(yk,order = c(p,0,q),include.mean = FALSE)  ## 'include.mean = TRUE', jos keskistystä ei tehdä (jolloin yk -> y1). 
fit1

## tai haluttu ARMA(p,q)x(P,Q)-kausivaihtelumalli 

fit1<-arima(yk,order = c(p,0,q),seasonal = list(order = c(P,0,Q), period = s),include.mean = FALSE) 
fit1

## Ennusteiden ja niiden luottamusvälien graafinen esittäminen

od <- options(digits = 4)  ## Valitaan (haluttaessa) ennusteiden keskivirheissä käytettyjen desimaalien määrä (tässä 4)
fcast <- predict(fit1, n.ahead = 16, se.fit = TRUE)  ## "n.ahead = 16" ilmoittaa laskettavieen ennusteiden lukumäärän (tässä 16, jolloin ennustetaan havintoperiodin loppuun), "se.fit = TRUE", kun keskivirheet palautetaan
fcast

f <- fcast$pred  ## Ennusteet vektoriin f
f <- f + mean   ## Lisätään ennusteisiin keskiarvo (sivuutetaan, jos keskistystä ei ole tehty)
se <- fcast$se  ## Ennusteiden keskivirheet vektoriin se

lv_ala <- f - 1.96*se  ## Vektori ennusteiden (tässä 95%:n) luottamusvälin alarajalle
lv_yla <- f + 1.96*se  ## Vektori ennusteiden (tässä 95%:n) luottamusvälin ylärajalle

forecast <- ts(c(y), start=c(1984,1), end=c(2012,4), freq=4) # Havaittu sarja ensin vektoriin forecast. 
plot(forecast, ylim=c(min(y, lv_ala, lv_yla),max(y, lv_ala, lv_yla)), lwd=1, lty=1,
   ylab="", main="US_GDP-sarja ja sen ennuste",xlab="")  ## Valitaan kuvan ylä- ja alaraja, jotta luottamusväli mahtuu kuvaan (voit joutua laajentamaan rajoja)

lines(f,lwd=2, col="blue", lty=5)  ## Piirretään ennuste sinisellä
lines(lv_ala, col="red", lwd=2, lty=4)  ## Piirretään ennusteen alaraja punaisella
lines(lv_yla, col="red", lwd=2, lty=4)  ## Piirretään ennusteen yläraja punaisella
abline(v=c(2009,1), lty=3, lwd=1)  ## Piirretään pystysuora viiva, joka erottaa havainnot ja ennusteet Ajankohdan määrittävän vektorin c(2009,1) valinnasta, ks. ?start


## Edellä "lwd" määrittelee viivan paksuuden, "lty" viivan tyypin ja "col" värin. Seuraavat vaihtoehdot ovat (ainakin) käytettävissä
## lwd: line width, positiivinen kokonaisluku oletusarvona 1
## lty: line type, vaihtoehtoina 0=blank, 1=solid (oletusarvo), 2=dashed, 3=dotted, 4=dotdash, 5=longdash, 6=twodash
## col: vaihtoehtoina "black", "red", "green3",  "blue" , "cyan", "magenta", "yellow", "gray"