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

## Luetaan mainittuun hakemistoon tallennettu aikasarja - tässä tapauksessa lake
## Nyt myös ensimmäisen havainnon alkamisajankohta (start=1892) ja havaintojen lkm/aikayksikkö (freq=1) on annettu

data=read.table("C:/Users/Erkka/Opiskelu/Tilastotiede/Stationaariset 2010/Laskuharjoitukset/lake.txt",header=TRUE)

Y=data
Y=ts(Y, start=1892,freq=1)

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

## 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) 
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
 
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ä tulostaman AIC:n arvo lasketaan eri tavalla kuin kurssilla ja alla.

n_parameters=length(fit$coef)
T=length(Yk)

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(Yk,order = c(p,0,q),include.mean = FALSE) ## 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
## viipymien määrä eli p+q (tai kausivaihtelumallissa p+q+P+q), jonka seuraavassa on oletettu olevan 2

resid=fit$residuals

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

## Vastaavasti neliöityjen residuaalien ja McLeod-Li -testin tapauksessa (nyt khi^2-jakauman vapausasteluvusta ei vähennetä viipymien määrää)

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")


## GARCH(r,s)-mallin parametrien estimointi ja diagnostiikka (s>0)
## Estimointi olettaa normaalijakauman eikä tulosta robusteja keskivirheitä

fit<-garch(dy,order = c(r,s)) ## estimoidaan GARCH(r,s,)-malli
summary(fit) ## estimointitulokset (a0=omega, ai=alpha_i, bi=beta_i)

## Jargue-Bera -testi testaa normaalisuusoletuksen realistisuutta

plot(fit) ## diagnostiikka graafisesti (näpäytä "Enter")