################
### PRVI CAS ###
################
data(AirPassengers)
AP <- AirPassengers
AP
class(AP) # odredjuje klasu objekta 
## Neki od metoda vezanih za klasu ts:
start(AP) # pocetak vs 
end(AP) # kraj vs 
frequency(AP) # period vs
## Genericke funkcije - mogu se primenjivati na objektima razlicitih klasa
summary(AP)
plot(AP,ylab="Passengers (1000s)")
## Komponente vremenske serije: 
  # trend - dugorocna pravilnost u kretanju vremenske serije, linearni ili polinom viseg reda, moze stohasticki (finansije)
  # sezonska komponenta - periodicna promena, period jedna godina ili uopsteno neki fiksni vremenski period
  # ponekad ciklusi koji ne odgovaraju nekom fiksnom periodu
# AP - razlozi za rastuci trend: bolji standard zivota nakon ww2, napredovanje avionske industrije, povecanje populacije
# AP - sezonska komponenta sa godisnjim periodom najvise putnika leti, najmanje krajem godine
# Na grafiku osim trenda/sezonske komonente nekad se mogu uociti i autlajeri/greske
# uklanjanje efekta sezonske varijacije, izravnavanje serije, tako da jasnije vidimo trend
plot(aggregate(AP)) # funkcijom aggregate sabira se po [frequency] elemenata, u nasem slucaju 12
boxplot(AP~cycle(AP))
#funkcija cycle odredjuje poziciju u ciklusu, 12 meseci-jedan ciklus
#sve vrednosti iz AP dele se u 12 kategorija i za svaku od njih se crta boxplot

#####
Maine <- read.csv("D:/.../Maine.txt", sep="")
attach(Maine)
class(Maine) 
# konvertovanje u ts:
Maine.month.ts <- ts(unemploy, start=c(1996,1), freq=12)
# sada pravimo sravnjenu seriju - od mesecne dobijamo godisnju
Maine.annual.ts <- aggregate(Maine.month.ts)/12
par(mfrow=c(1,2))
plot(Maine.month.ts)
plot(Maine.annual.ts)
par(mfrow=c(1,1))
(Maine.Feb <- window(Maine.month.ts, start = c(1996,2), freq = TRUE))
(Maine.Aug <- window(Maine.month.ts, start = c(1996,8), freq = TRUE))
# funkcija window - primer
# a<-1:100
# window(a,start=1, deltat=5)
# ako se zadaje vs vraca neki podskup te vs, moze se zadati i end
(Feb.ratio <- mean(Maine.Feb) / mean(Maine.month.ts))
(Aug.ratio <- mean(Maine.Aug) / mean(Maine.month.ts))
# nezaposlenost veca u februaru za 22% od proseka, a u avgustu 18% manja od proseka
US.month <- read.csv("D:/.../USunemp.txt", sep="")
attach(US.month)
US.month.ts <- ts(US.month, start=c(1996,1), freq = 12)
plot(US.month.ts, ylab = "unemployed (%)")
# odavde vidimo slicnosti izmedju serija - opadanje nezaposlenosti pocetkom '00
# razlika je da u Maine-u ne opada nezaposlenost oko 2006. kao u ostatku SAD

#####
CBE <- read.delim("D:/.../cbe.txt")
CBE[1:4, ] # prve 4 vrste i sve kolone
class(CBE)
#ne mora da se navede drugi argument za start ili end, onda se pretpostavljaju 1 odnosno 12
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts <- ts(CBE[, 1], start = 1958, freq = 12)
# sa funkcijom cbind unutar plot, dobijamo vise vremenskih serija
plot(cbind(Elec.ts, Beer.ts, Choc.ts))
# primecujemo rastuci trend za sva tri, posledica povecanja populacije
# ovo je multiple vremenska serija
AP.elec <- ts.intersect(AP, Elec.ts) # dve serije spajamo u jednu dvodimenzionu na preseku vremena
start(AP.elec)
end(AP.elec)
AP.elec[1:3,]
AP <- AP.elec[,1] 
Elec <- AP.elec[,2] 
par(mfrow=c(1,2))
plot(AP, main = "", ylab = "Air passengers / 1000's")
plot(Elec, main = "", ylab = "Electricity production / MkWh")
# slicnost izmedju ovih vs i po rastucem trendu i po sez kom (struja se proizvodi zimi vise)
layout(1:1)
plot(as.vector(AP), as.vector(Elec),
     xlab = "Air passengers / 1000's",
     ylab = "Electricity production / MWh")
#potrebno je pretvoriti serije u nizove
abline(reg = lm(Elec ~ AP))
#vidimo da su linearno zavisne, preko korelacije (koja meri lin zav) dobija se
cor(AP,Elec) # cor genericka, prima i klase ts
# medjutim iako su korelisane ne znaci da je jedna posledica druge
# AP u USA, Elec u AUS

#####
# exchange rate - koliko jedne valute se dobija za jedinicnu vrednost druge (r), ili obrnuto (bice 1/r)
Z <- pounds_nz <- read.csv("D:/.../pounds_nz.txt", sep="")
Z.ts <- ts(Z,frequency=4,start=1991)
plot(Z.ts, xlab="time (years)",ylab="quarterly exchange rate in $NZ/pound")
# Mozemo izdvojiti kratkorocne (short-term) trendove, npr pad od 1992 do 1996
# pa rast od 1996 do 1998, ali smer kretanja trenda (opadajuci, rastuci) se menja
# u nepredvidivim trenucima. Ovakvi trendovi se zovu STOHASTICKI trendovi
Z.92.96<-window(Z.ts,start=1992,end=c(1996,1))
Z.96.98<-window(Z.ts,start=1996,end=c(1998,1))
#nismo naveli frequency=TRUE, pa uzima sve vrednosti u navedenom intervalu
par(mfrow=c(1,2))
plot(Z.92.96)
plot(Z.96.98)
par(mfrow=c(1,1))
# Da smo hteli na pocetku 1996. da predvidimo dalje kretanje
# kursa,bez nekih dodatnih informacija, samo na osnovu podataka,
# napraviili bismo gresku. Predvidjanje je problem kod stohastickih trendova.
# Postoje testovi kojima se utvrdjuje da li je trend stohasticki.

#####
global <- scan("D:/.../global.txt",  sep="") #ucita u niz
Global.ts <- ts(global,frequency=12,start=c(1856,1),end=c(2005,12))
start(Global.ts)
end(Global.ts)
# najvise nas zanima trend pa izravnjujemo seriju
# moguce je navesti funkciju koju cemo primeniti na skup spojenih podataka, a da nije samo sabiranje
Global.annual <- aggregate(Global.ts, FUN=mean)
plot(Global.ts)
plot(Global.annual)
#izdvajamo samo onaj period kritican za globalno zagrevanje
New.series <- window(Global.ts, start=c(1970,1))
plot(New.series)
New.time <- time(New.series)
# funkcija time u niz stavlja sve vremenske trenutke u kojima su vrsena merenja
# godina+(mesec-1)/12
# malopre smo imali regresiju kojom smo ispitivali odnos dva obelezja
# sada trend ocenjujemo regresionom pravom
# i sto se tice predvidjanja imamo opravdanja jer ovakav trend, zbog globalnog zagrevanja nije slucajan
abline(lm(New.series~New.time))

### MODELI ###

## Aditivni model ##
# x_t= m_t + s_t + z_t
# m_t - trend
# s_t - 
# z_t - proces greske, u opstem slucaju korelisane slucajne velicine sa srednjom vrednostu 0
## Multiplikativni model ##
# x_t=m_t*s_t+z_t
# ako je sezonska komponenta raste kako trend raste 
# x_t=m_t*s_t*z_t 
# Svodi se na aditivni model kada logaritmujemo (serija mora da bude pozitivna)
# problem kada korsitimo funkciju exp() da bismo se vratili na originalni model

# Kako oceniti trend? 
# do sada smo pokusali regresijom, funkciju aggregate smo koristili da bismo iskljucujuci pretpostavljenu godisnju sezonsku komponentu
# dobili bolju sliku trenda

## Metod pokretnih proseka ##
# metod izravnanja serije, za ocenu trenda u tacki t
# koristi se prosek nekoliko elemenata u okolini te tacke (za svako t, pokretno)
# Ako pretpostavimo da je sezona jedna godina, da bi se anulirala
# sezonska komponenta treba sabrati svih 12 meseci, ali ako ocenjujemo trend u julu
# trebalo bi da se odlucimo da li cemo sabrati 5 meseci pre jula i 6 posle
# ili 6 pre i 5 posle. Da bismo izbegli zabunu resenje je sredina izracunatih
# vrednosti u oba slucaja 
# odavde se dobije ocena sezonske komponente oduzimanjem odnosno deljenjem pocetne serije

# Dekompozicija
# za ocenu trenda koristi se metod pokretnih proseka 
Elec.decom.add <- decompose(Elec.ts)
plot(Elec.decom.add)
Elec.decom.mult <- decompose(Elec.ts, type = "mult")
plot(Elec.decom.mult)
# na osnovu prvog grafika vidimo da je pogodniji multiplikativni model jer sezonska komponenta raste sa porastom trenda
# kad se u funkciji decompose navede argument type="mult" misli se na multiplikativni x_t=m_t*s_t*z_t model
Trend <- Elec.decom.add$trend
Trend #Vidimo da nedostaje prvih 6 i poslednjih 6 elemenata
Seasonal <- Elec.decom.add$seasonal
ts.plot(cbind(Trend,Trend+Seasonal),lty=2)
#####
# Programiramo funkciju za ocenjivanje trenda metodom pokretnih proseka
# rezultat treba da bude isti kao trend dobijen funkcijom decompose
MA<-function(ts){
  x<-vector()
  ts1<-as.vector(ts)
  x[1:6]<-NA; x[(length(ts1)-6):(length(ts1))]<-NA
  for(i in 7:(length(ts1)-6)) 
  {
    x[i]<-(ts1[i-6]/2+ts1[i+6]/2+sum(ts1[(i-5):(i+5)]))/12
  }
  x.ts<-ts(x,start=c(start(ts)[1],start(ts)[2]),freq=12)
  return(x.ts)
}
plot(Trend)
plot(MA(Elec.ts))
# Jos je potrebno videti kako se tacno dobija sezonska komponenta.
# Za detrendovanu vremensku seriju izracunaju se srednje vrednosti po svakom mesecu, pa se taj niz, od 12 vrednosti centrira
SV<-function(ts,trend){
  m<-t(matrix(data = ts-trend, nrow = 12))
  seasonal<-colMeans(m, na.rm = TRUE)
  seasonal<-seasonal-mean(seasonal)
  ts.seasonal<-ts(rep(seasonal,end(ts)[1]-start(ts)[1]+1),start=start(ts)[1],freq=12)
  plot(ts.seasonal,ylab="Seasonal Variation")
  return(seasonal)
}
SV(Elec.ts,Trend)
plot(Seasonal)
# na osnovu ova dva grafika (ista su ) uverili smo se kako tacno radi
# funkcija decompose za aditivni model, a slicno je i za multiplikativni
# Rucno smo dobili trend i sezonsku varijaciju, a gresku dobijamo kada
# od originalne serije oduzmemo trend i sezonsku komponentu