####################
### ARIMA MODELI ###
####################

### NESEZONSKI ARIMA MODELI ###

## PRIMER1 ##
#install.packages("fpp")
#library(fpp)
#install.packages("forecast")
#library(forecast)
plot(usconsumption[,1]) # Promene u procentima, potrosnje u SAD, kvartalna serija
abline(v=1970:2010,lty="dotted")
#***zakljucivanje sa grafika***
# ne izgleda kao da ima sezonsku komponentu
# niti kao da je potrebno diferenciranje, nema ociglednog trenda
usc <- usconsumption[,1]
tsdisplay(usc)
#***transformisanje***
# necemo logaritmovati, niti transformisati, vec se posmatraju
# promene samo u procentima, sto je zgodno da bude stacionaran proces
#***prvi model***
# na osnovu ACF i PACF probacemo MA(3) model, a mozemo i AR(3)
fit1 <- Arima(usc, order=c(3,0,0))
fit2 <- Arima(usc, order=c(0,0,3))
fit1$aic
fit2$aic # malo je bolji AR(3) model
#***automatsko odredjivanje reda***
fit <- auto.arima(usconsumption[,1],seasonal=FALSE)
fit
#***popravljanje***
tsdisplay(residuals(fit1))
tsdisplay(residuals(fit2)) # i reziduali slicno izgledaju
#***predvidjanje***
# odaberemo ARIMA(3,0,0)
plot(forecast(fit1,h=10)) # prikazani intervali poverenja su 95%-ni i 80%-ni

## PRIMER2 ##
# paket fpp
class(oil) # godisnja proizvodnja nafte u Saudijskoj Arabiji
plot(oil)
adf.test(oil) # nije stacionaran
kpss.test(oil)
ndiffs(oil) # ocekivano vraca 1
doil <- diff(oil) # jednom diferenciramo
tsdisplay(doil)
fit1 <- Arima(oil, order=c(0,1,0))
fit1$aic
res <- residuals(fit1)
tsdisplay(res) # nema znacajnih autokorelacija
auto.arima(oil) # i ovim smo dosli do istog modela
plot(forecast(fit1,h=6))

## PRIMER3 ##
# paket astsa
class(oil) # cene WTI - West Texas Intermediate (tip sirove nafte) cija cena dalje utice i na cene preradjevina
#cena je u dolarima po barelu serija je nedeljna od 2000. do sredine 2010. god
frequency(oil) # 52 jer je nedeljna vs
plot(oil)
abline(v=2000:2010,lty="dotted") # ne vidimo sezonsku komponentu
# serija ne izgleda stacionarno, i izgleda da disperzija raste sa porastom trenda
# prvo logaritmujemo
logoil <- log(oil)
plot(logoil) # smirili smo disperziju
adf.test(logoil)
kpss.test(logoil) # napomena: ne moze se uvek verovati testovima jedinicnog korena, ovde ipak mozemo
ndiffs(oil) # ocigledno je da serija nije stacionarna, pa cemo diferencirati jednom
dlogoil <- diff(logoil)
tsdisplay(dlogoil)
# zbog velikog broja autokorelacija, nije jasno sta se desava tacno
# acf i pacf izgledaju slicno
# imamo pozitivnu acf(1) pokusacemo dodavanjem jednog ar parametra da ubije
# pozitivnu autokorelaciju
fit1 <- Arima(oil,order=c(1,1,0),lambda=0)
res <- residuals(fit1)
tsdisplay(res) 
# sada cemo dodati jedan MA parametar jer imamo negativnu prvu autokorelaciju
fit2 <- Arima(oil,order=c(1,1,1),lambda=0)
res <- residuals(fit2)
tsdisplay(res)
#ovim smo prvih nekoliko acf i pacf vrednosti smirili, tako da su u granicama
Box.test(res,type="Ljung-Box",lag=104, fitdf=2)
# prolazi test! iako ima i znacajnih autokorelacija
# sveukuono ih ima puno, pa je ok da neke ne budu u intervalu
qqnorm(res) # i ne izgleda bas normalno ali nema veze
auto.arima(oil,d=1,lambda = 0,seasonal=FALSE) #dobijamo isti model
# predvidjanje
plot(forecast(fit2,h=52))
# skoro da je ravna linija, predvidjanje za narednu godinu
# nema puno smisla jer vidimo kako se menja puno i tokom nedelja
# ipak znace nam intervali poverenja, kao donja i gornja granica
# intervala u kom ce se kretati serija

### SEZONSKI ARIMA MODELI ###

# primer ACF i PACF funkcija nekog sezonskog ARIMA modela
phi=c(rep(0,11),.8)
ACF=ARMAacf(ar=phi, ma=-0.5, 50)[-1]
PACF=ARMAacf(ar=phi, ma=-0.5, 50, pacf=TRUE)
plot(ACF, type="h")
abline(h=0)
plot(PACF, type="h")
abline(h=0)

# PRIMER1 #
class(euretail) # Quarterly retail trade index in the Euro area (17 countries), 1996-2011
frequency(euretail)
#***zakljucivanje sa grafika***
plot(euretail)
abline(v=1996:2011,lty="dotted") # Vidi se i sezonska komponenta
#***diferenciranje***
nsdiffs(euretail) # prvo sezonski diferenciramo
tsdisplay(diff(euretail,4))
ndiffs(diff(euretail,4)) # ni ovo ne izgleda stacionarno, pa cemo i obicno diferencirati
tsdisplay(diff(diff(euretail,4))) # ovo izgleda malo stacionarnije, srednja vrednost prblizno nula
# na ACF i PACF vide se znacajnost na sezonskim zadrskama, sto nam vise i odgovara
# odabrali smo d=1,D=1
#***odabir ostalih parametara***
# imamo negativnu vrednost ACF(4) i pozitivnu ACF(1)
#P pobacemo sa ARIMA(0,1,1)(0,1,1)_4 modelom, pa mozda budemo i dodali jos param
# funkcija Arima, slicno kao arima, ima jos neke mogucnosti koje necemo koristiti
# arima je iz stats paketa, Arima iz forecast
fit <- Arima(euretail, order=c(0,1,1), seasonal=c(0,1,1))
tsdisplay(residuals(fit))
# probati i 1,1,0 0,1,1
#***popravljanje***
# Po rezidualima vidimo da smo odabranim SMA=1 resili sezonske parametre
# ali znacajne korelacije na manjim zadrskama sugerisu da treba dodati jos neki nesezonski
# Pitanje je koji parametar dodati. U odeljku Odabir novog modela, kod nesezonskih arima
# receno je da ako je (negativni) siljak na ACF za male lag, onda se dodaje jedan MA par
# a ako je (pozitivni) siljak na PACF za male lag onda se dodaje jedan AR par
# Prema tome ovde bismo mogli da dodamo i jedan AR nesezonski, ili MA nesezonski
# (dodavanjem MA jos dobijamo pomenuti uobicajeni sezonski model)
# Nema jedinstvenog pravila, mozete probati i dodavanjem AR parametara
# Mi cemo dodati MA
fit4 <- Arima(euretail, order=c(1,1,1), seasonal=c(0,1,1))
res <- residuals(fit4)
tsdisplay(res)
fit4$aic
fit2 <- Arima(euretail, order=c(0,1,2), seasonal=c(0,1,1))
res <- residuals(fit2)
tsdisplay(res)
fit2$aic
Box.test(res, lag=16, fitdf=3, type="Ljung") # lag=16, znaci posmatramo ukupno 16 autokor
# fitdf=3 jer smo ocenili 3 parmetara 2xAR, 1XSMA
# nulta hipoteza je da su autokorelacije 0
qqnorm(res)
# Ovo je ok model gledajuci po p vrednosti Ljung Box testa
# ali ipak vidimo jos znacajne ACF 3 i PACF 3
# probacemo dodavanjem jos jednog MA parametra
fit3 <- Arima(euretail, order=c(0,1,3), seasonal=c(0,1,1))
res <- residuals(fit3)
tsdisplay(res)
Box.test(res, lag=16, fitdf=4, type="Ljung")
fit3$aic
# Sada i reziduali bolje izgledaju, i veca je p vrednost testa, a mozemo da uporedimo
# AIC vrednosti, ili AICc
fit2$aicc
fit3$aic #ovo je bolje jer je manje
#***Predvidjanje***
# Nastavicemo da radimo sa paketom forecast, u kojoj je funkcija za predvidjanje forecast
# inace postoji i predict funkcija
plot(forecast(fit3,h=12)) # predvidjanje za naredne 3 godina
# Predvidjanja prate skorasnji opadajuci trend, ali intervali poverenja su veliki
# i dozvoljavaju da serija u nekom trenutku postane i rastuca
#***Automatsko odredjivanje modela***
# Pomenuta funkcija auto.arima moze odrediti i sezonske redove
auto.arima(euretail)
# Nakon malo cekanja, dobije se da je najbolji model ARIMA(1,1,2)(0,1,1)
# a njegova vrednost AICc je cak veca od naseg odabranog modela
# To se desava jer ova funkcija ne proverava bas sve modele, nego pravi neke precice
# u nalazenju modela, pa nekad i ne da stvarno najbolji model
# Te precice se mogu iskljuciti, pa sledecm pozivom dobijamo najbolji model
# koji smo i sami identifikovali
auto.arima(euretail,stepwise=FALSE, approximation=FALSE)
# Vidimo da serija nije stacionarna

# PRIMER2 #
#***zakljucivanje sa grafika***
# Mesecna serija, prodaja lekova sa kortikosteriodima u Australiji
plot(h02)
# sezonska komponenta je jasna i izgleda kao da raste sa porastom trenda
# serija je pozitivna pa je mozemo logaritmovati
lh02 <- log(h02)
par(mfrow=c(2,1))
plot(h02, ylab="H02 sales (million scripts)", xlab="Year")
plot(lh02, ylab="Log H02 sales", xlab="Year")
# Sa ovog grafika se moze videti kako smo logaritamskom transformacijom
# malo smirili disperziju (oko 2005 godine)
# i dalje imamo blago rastuci trend i ociglednu sezonsku komponentu
#***diferenciranje***
tsdisplay(diff(lh02,12), 
          main="Seasonally differenced H02 scripts", xlab="Year") # prvo sezonsko
# ovo nam je mogla preporuciti i funkcija
nsdiffs(lh02)
# da li je potrebno i obicno diferencirati?
kpss.test(diff(lh02,12)) # kpss test kaze da, a kod adf je p vrednost 0.06
# to takodje kaze i funkcija ndiffs
ndiffs(diff(lh02,12)) 
# pogledajmo grafik ACF i PACF kada i obicno diferenciramo
tsdisplay(diff(diff(lh02,12),1))
# receno je (vidi pdf odeljak Prvi korak u identifikaciji modela) da ako su ACF
# vrednosti manje od -0.5 onda je previse puta diferencirano
# Zbog tog ipak necemo jos jednom diferencirati
tsdisplay(diff(lh02,12))
# ovo odgovara drugom opisanom uobicajenom sezonskom modelu ARIMA(p,0,0)(0,1,1)+c
# imamo dosta pozitivnih autokorelacije koje ubijamo dodavanjem AR parametara
fit1<-Arima(h02, order=c(1,0,0),seasonal=c(0,1,1),include.mean=TRUE,lambda=0)
res<-residuals(fit1)
tsdisplay(res)
# od mnogo pozitivnih, dodavanjem AR parametara, dobili jako male ACF i PACF 1
#***popravljanje***
# moracemo da popravimo ovu situaciju,ali ne dodavanjem AR parametara jer cemo tako
# napraviti samo jos negativnije ACF, mozete to i probati
# dakle, odustajemo od uobicajenog modela
# probacemo dodavanjem MA nesezonskog parametara 
fit2<-Arima(h02, order=c(1,0,1),seasonal=c(0,1,1),lambda=0)
res<-residuals(fit2)
tsdisplay(res)
fit3<-Arima(h02, order=c(1,0,2),seasonal=c(0,1,1),lambda=0)
res<-residuals(fit3)
tsdisplay(res) # nikako se ne popravlja situacija, na velikim lag znacajne autokorelacije
# promenicemo i sezonske parametre
# krecemo otpocetka
tsdisplay(diff(lh02,12))
# probacemo sa dodavanjem jos jednog SMA parametra
# stavicemo jedan AR parametar zbog pozitivnih autokorelacija na malim lag
fit4<-Arima(h02, order=c(1,0,0),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit4)
tsdisplay(res)
fit4$aicc
# dodajemo jedan MA zbog sada negativnih autokorelacija na malim lag
fit5<-Arima(h02, order=c(1,0,1),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit5)
tsdisplay(res)
fit5$aicc
# dodajemo jos dva AR cime se popravlja situacija na malim lag
# iako ostaju neke znacajne ACF i PACF za velike lag
# ali to nije na sezonskim zadrskama pa ignorisemo
fit6<-Arima(h02, order=c(3,0,1),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit6)
tsdisplay(res)
fit6$aicc
Box.test(res, lag=36, fitdf=7, type="Ljung")
# uzimamo za lag (broj prvih autokorelacije koje se zajedno posmatraju) da je neki
# umnozak sezonskog perioda, fitdf je broj ocenjenih parametara, a to je 7 zbog mean
# Reziduali naseg najboljeg modela nisu prosli test
# Nekad nije moguce naci model koji ce proci sve testove
# probacemo sa auto.arima
#***automatsko odredjivanje reda***
# nametnucemo nase diferenciranje i transformaciju
fit<-auto.arima(h02,lambda=0,d=0,D=1,max.order=9,stepwise=FALSE, approximation = FALSE)
# max.order je najveca vrednost za p+q+P+Q
# ova funkcija se duze ceka, dobija se da je najbolji red ARIMA(4,0,3)(0,1,1)
# i ovaj model prolazi Ljung Box test
Box.test(residuals(fit),lag=36,fitdf=8,type=Ljung)
#***moc predvidjanja kao mera kvaliteta modela***
#Funkcija getrmse radi sledece
getrmse <- function(x,h,...){
  train.end <- time(x)[length(x)-h]
  test.start <- time(x)[length(x)-h+1]
  train <- window(x,end=train.end)
  test <- window(x,start=test.start)
  fit <- Arima(train,...)
  fc <- forecast(fit,h=h)
  return(accuracy(fc,test)[2,"RMSE"])
}
# za odabrano h, izdvaja poslednjih h vrednosti iz serije
# to je trening skup za nas model, na tom skupu za dati red modela se ocenjuju koef
# test skup su poslednjih h vrednosti serije, na tom skupu testiramo nas model
# tako sto racunamo koja je greska predvidjanja
# Mi imamo stvarne vrednosti i poredimo ih sa predvidjenima
# Povratna vrednost funkcije je koren srednje kvadratne greske (koren sume kvadrata greski)
# probamo za razne modele da vidimo kako se krecu vrednosti
getrmse(h02,h=24,order=c(3,0,0),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,2),seasonal=c(2,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(1,1,0),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(0,1,2),lambda=0) # nas najbolji
getrmse(h02,h=24,order=c(3,0,1),seasonal=c(1,1,1),lambda=0)
getrmse(h02,h=24,order=c(4,0,3),seasonal=c(0,1,1),lambda=0) # najbolji prema auto.arima
getrmse(h02,h=24,order=c(3,0,3),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(4,0,2),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(3,0,2),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,3),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,4),seasonal=c(0,1,1),lambda=0)
getrmse(h02,h=24,order=c(2,1,5),seasonal=c(0,1,1),lambda=0)
# ne menjaju se puno, uglavnom onim modelima sa manjim AIC odgovara i bolje predvidjanje
# ali nas model je po predvidjanju bolji od auto.arima modela
# a auto.arima model prolazi test reziduala, dok nas nije
# Bitno je da kada se porede modeli prema vrednosti AIC oni moraju biti diferencirani
# isti broj puta
# Koriscenje test i trening skupa je drugi nacin da se porede modeli
# Tu nije bitno kog je reda diferenciranje
#***predvidjanje***
# Definitivno smo odabrali ARIMA(3,0,1)(0,1,2) model
# Predvidjamo vrednosti 
plot(forecast(fit6)) 
# po defaultu predvidjanje je 2*sezonski period, dakle 2 godine