#################
### TRECI CAS ###
#################

### HOLT-WINTERS METOD 

# Metod koji koristimo kada imamo i trend i sezonsku komponentu, moze da se koristi i za predvidjane
# Jednacine ADITIVNOG HOLT-WINTERS modela
wine <- read.delim("D:/.../wine.txt")
attach (wine)
# prodaja australijskog vina (u 1000l) od 1.1980. do 7.1995. za razlicite vrsta vina (multiple ts)
sweetw.ts <- ts(sweetw, start = c(1980,1), freq = 12)
plot(sweetw.ts, xlab= "Time (months)", ylab = "sales (1000 litres)")
# prodaja slatkog belog vina belezi veliki rast nakon 1985. god.
sweetw.hw <- HoltWinters(sweetw.ts)
sweetw.hw # ispisuju se parametri modela i koeficijenti (vrednosti koje sluze za predvidjanje)
sweetw.hw$SSE # minimizirana greska predvidjanja jedan korak unapred
head(sweetw.hw$fitted) # pocinje sa godinom zakasnjenja, moguce je ocenjivat nakon s (ovde 1 godina)
plot(sweetw.hw) # crna - originalna serija; crvena - ocenjena serija iz prethodnog reda
plot (sweetw.hw$fitted) # crvena linija na prethodnom grafiku je xhat u ovom grafiku
# xhat je zbir level+trend+seasonal
sweetw.hw$coef # niz ocenjenih koeficijenata modela u poslednjim trenucima
# a=a_n, b=b_n, s_{n-p+1},...,s_n su ocena 12 sezonskih komponenti
# Na osnovu njih se vrsi predvidjanje 
# Jednacina predvidjanja za k koraka

# Jednacine MULTIPLIKATIVNOG HOLT-WINTERS modela
sweetw.hw <- HoltWinters(sweetw.ts, seasonal="mult")
sweetw.hw
plot(sweetw.hw)
sweetw.hw$SSE
# Vidimo da je beta=0, to znaci da se za vrednost nagiba uzima poslednja ocenjena vrednost
# Znaci, ocenjeno je da se nagib trenda sporo menja, za razliku od nivoa i sezonske komponente
# multiplikativni model je bolji od aditivnog - PROVERITI, za njega cemo raditi predvidjanje
# Osnovna funkcija za prognoziranje - predict
sweetw.hw.predict <- predict(sweetw.hw, n.ahead=4*12)
ts.plot(sweetw.ts,sweetw.hw.predict,lty=1:2)
# Drugi nacin za prognoziranje- paket forecast
library(forecast)
sweetw.hw2 <- forecast(sweetw.hw, h=4*12)
# argument h - za koliko vrednosti unapred predvidjamo
plot(sweetw.hw2)
# Svejedno je da li se koristi funkcija predict ili forecast - dobijaju se iste prognozirane vrednosti
# Samo kod funkcije forecast dobijamo i intervale poverenja (kako??)
# U ovom primeru nismo se bavili rezidualima - tj ispitivanjem da li je model dobar
# Videti u narednim primerima sta se proverava za reziduale pa testirajte to i na ovom primeru


### EKSPONENCIJALNO IZRAVNANJE

Motor <- read.csv("D:/.../motororg.txt", sep="")
attach(Motor)
# Broj zalbi upucenih AMS-u u periodi od 1996. do 1999. godine
# Zelimo da predvidimo broj zalbi u 2000. i ocenimo trend i sezonsku komponentu
Comp.ts <- ts(complaints, start = c(1996, 1), freq = 12)
plot(Comp.ts, xlab = "Time / months", ylab = "Complaints")
# na osnovu grafika vidimo da je srednja vrednost priblizno 20,
# ali ne vidimo ni jasno trend ni sezonsku varijaciju,
# zato cemo koristiti eksponencijalno ravnavanje
# Potrebno je podesiti ove parametre beta i gamma na FALSE
# a posto se ne navodi vrednost argumenta alpha, ona se ocenjuje
# tako da suma kvadrata gresaka predvidjanja za jedan korak bude minimalna
Comp.hw1 <- HoltWinters(complaints, beta = FALSE, gamma = FALSE) 
Comp.hw1
plot(Comp.hw1)
# Ako sami zadamo alpha:
# Za alpha=1 - ista serija, pomerena jedan korak unapred
# jer nemamo a_{1-1} (a_0) za ocenjivanje a_{1}, pa izravnjena serija kod eksp izravnavanja
# krece od a_2
# Kada imamo i sezonsku komponentu, godisnju npr, serija je pomerena celu godinu unapred
Comp.hw2 <- HoltWinters(complaints, alpha = 0.3, beta=FALSE, gamma=FALSE)
plot(Comp.hw2)
Comp.hw1$fitted # Fitovane vrednosti, tj. vrednosti izravnate serije su bas jednake nivou
(Comp.hw1.predict <- predict(Comp.hw1, n.ahead=12))
# 17.7 je ocena na kraju 1999 godine (za bilo koji korak unapred)

rain <- read.delim("D:/.../rain.txt")
attach(rain)
# U seriji su date godisnje vrednosti kolicine kise (u incima) u Londonu od 1813
rainseries <- ts(Rainfall,start=c(1813)) 
plot(rainseries)
# Ne vidimo jasan trend ni neki ciklus koji se ponavlja, nivo je priblizno konstatan
rainseriesforecasts <- HoltWinters(rainseries, beta=FALSE, gamma=FALSE)
rainseriesforecasts # Dobija se jako mala vrednost alpha
rainseriesforecasts$fitted
# Uticaj razlicitih vrednosti za alpha kod eksponencijalnog izravnavanja
plot(rainseriesforecasts) # Nacrtana pocetna, i izravnata serija
plot(HoltWinters(rainseries, alpha=0.25, beta=FALSE, gamma=FALSE), main="Alpha=0.25")
plot(HoltWinters(rainseries, alpha=0.5, beta=FALSE, gamma=FALSE), main="Alpha=0.5")
plot(HoltWinters(rainseries, alpha=0.75, beta=FALSE, gamma=FALSE), main="Alpha=0.75")
plot(HoltWinters(rainseries, alpha=1, beta=FALSE, gamma=FALSE), main="Alpha=1")
library(forecast)
(rainseriesforecasts2 <- forecast(rainseriesforecasts, h=50))
plot(rainseriesforecasts2)
# vidimo da su sve prognoze iste, jednake poslednjoj vrednosti nivoa
# ovde su date i granice intervala poverenja za ove ocenjene vrednosti
# granice intervala poverenja racunaju se pod pretpostavkom da su greske nezavisne,
# normalno raspodeljene sa konstantnom disperzijom 

# Kako rucno dobiti 95% interval poverenja?

# Pozabavimo se sada serijom reziduala u cilju provere da li je HW odgovarajuci model
# na osnovu koga predvidjamo. 
# Treba proveriti sledece:

# 1. Da li nema serijske korelacije u seriji reziduala-korelogram?

acf(rainseriesforecasts2$residuals[-1], lag.max=20)
# Testiranje znacajnosti korelacija moze se sprovesti, osim testom datim na grafiku korelograma
# jos jednim testom. To je Ljung-Box test koji takodje testira da li ima neke statisticki znacajne korelacije
Box.test(rainseriesforecasts2$residuals, lag=20, type="Ljung-Box")
# Nulta hipoteza je da su sve korelacije jednake 0, i posto je ovde p-vrednost testa velika
# prihvatamo nultu hipotezu

# 2. Da li su reziduali normalno raspodeljeni?

qqnorm(rainseriesforecasts2$residuals[-1])
# normalnost reziduala mozemo proveriti na vise nacina, na osnovu qqplota, gde poredimo
# kvantile serije reziduala i teorijske kvantile normalne raspodele 

#3. Da li je disperzija serije reziduala konstantna tokom vremena?

# Moze se zakljuciti na osnovu grafika
# Sa grafika bismo rekli da je mozda ipak manja disperzija na pocetku, a veca kasnije
plot(rainseriesforecasts2$residuals[-1])

skirts <- read.delim("D:/.../skirts.txt")
attach(skirts) # Godisnja serija precnik suknji, od 1866 god
skirts.ts <- ts(Diameter,start=c(1866))
plot(skirts.ts)
# sa grafika ne vidimo sezonsku komponentu, ali postoje neki kratkorocni trendovi
# na pocetku rastuci posle opadajuci
# Primenicemo HW model, samo bez parametra \gamma i sezonske komponente
skirts.ts.hw <- HoltWinters(skirts.ts, gamma=FALSE)
skirts.ts.hw
# vidimo da je beta=1, dakle nagib se ocenjuje kao razlika vrednosti serije u uzastopnim trenucima
# rekli smo da je to dobra ocena kad je trend nepredvidiv
library(forecast)
skirts.ts.hw.forecast <- forecast(skirts.ts.hw, h=19)
skirts.ts.hw.forecast
plot(skirts.ts.hw.forecast)
skirts.ts.hw.predict <- predict(skirts.ts.hw, n.ahead= 19) 
ts.plot(skirts.ts, skirts.ts.hw.predict, lty = 1:2) 
acf(skirts.ts.hw.forecast$residuals[-(1:2)],lag.max=20)
# kao argument acf mozemo zadati koliko vrednosti autokorelacione funkcije zelimo sa lag.max 
# provera da li je model odgovarajuci na osnovu korelograma: imamo jednu stat znacajnu korelaciju
# sto je ok, jer treba je dozvoljeno da 5% korelacija upadne u kriticnu oblast (tj 1 korelacija ako gledamo njih 20)
# Provera Ljung-Box testom potvrdjuje da je model dobar
Box.test(skirts.ts.hw.forecast$residuals, lag=20, type="Ljung-Box")
plot(skirts.ts.hw.forecast$residuals[-(1:2)])
#Provera da li je disperzija konstantna tokom vremena - na osnovu grafika rekli bismo da jeste

data(AirPassengers) # serija analizirana na prvom casu
# Posto se cinilo da se sezonski efekat povecava kako trend raste, koristicemo multiplikativni HW model
AP <- AirPassengers
AP.hw <- HoltWinters(AP, seasonal = "mult")
AP.hw
plot(AP.hw)
AP.predict <- predict(AP.hw, n.ahead= 4 * 12)
ts.plot(AP, AP.predict, lty = 1:2) 
plot(AP.predict)