#################
### DRUGI CAS ###
#################
# Beli sum - Svaki proces sa nezavisnim vrednostima, srednje vrednosti 0 i disperzije sigma^2
# Nije obavezno da bude iz normalne raspodele, ako jeste iz normalne raspodele onda se to zove Gausov beli sum
# Ovo je idealni model za reziduale
# Nakon sto uklopimo seriju u neki model (modele cemo kasnije definisati) idealno je da razlika izmedju realizovane i modelom predvidjene vrednosti bude beli sum
w <- rnorm(100)
plot(w, type = "l")
x <- seq(-3,3, length = 1000)
hist(rnorm(100), prob = TRUE)
points(x, dnorm(x), type = "l") # funkcijom points na vec postojeci grafik dodajemo tacke
#ovo je drugi nacin za dodavanje linije na grafik
curve(dnorm(x),xlim=c(-3,3),add=TRUE) 
w <- rnorm(500,0,1) # 500 N(0,1) realizacija
# MA belog suma - izravnavanje serije
v <- filter(w, sides=2, rep(1/3,3)) # ova fja je mogla umesto rucno napravljene funkcije MA na prvom casu
par(mfrow=c(2,1))
plot.ts(w, main="white noise")
plot.ts(v, main="moving average") #vidimo da je druga serija mirnija po ponasanju

# ako nekad dobijate gresku: Error in plot.new() : figure margins too large
# par(mar=c(1,1,1,1))

# Braunovo kretanje (Vinerov proces)
# U finansijskoj matematici jedan moguci model za cene akcija je proces exp(W_t) gde je W_t Braunovo kretanje
# ovo je geometrijsko Braunovo kretanje
# za simuliranje koristi se B_0=0, svojstvo nezavisnosti prirastaja i raspodela B_t-B_s
# ako hocemo da simuliramo ograniceni smo diskretnim vremenom
# u 1000 tacaka aproksimirati Braunovo kretanje na intervalu [0,1000]
# prirastaji na vremenskom intervalu duzine k imaju N(0,k) raspodelu
x <- w <- rnorm(1000)
for (t in 2:1000) x[t] <- x[t - 1] + w[t]
plot(x, type = "l")
# Braunovo kretanje na [0,1]
# posto simuliramo na vremenskom intervalu [0,1] interval izmedju dve uzastopne
# opservacije uzimamo da bude mali, npr 0.001, onda je to i disperzija razlike te dve opservacije
t <- seq(0,1,0.001) #1001 podeona tacka
B <- vector()
B[1] <- 0
for(i in 1:1000){
  B[i+1] <- B[i]+rnorm(1,0,0.001)
}
plot(t,B, type="l")
# sada prirastaj izmedju susednih tacaka ima raspodelu N(0,0.001), a malopre je imao N(0,1)
# dakle nisu samo skalirane vrednosti na x-osi, a to se moze videti i po vrednostima na y osi, kod slucajnog lutanja su manje, jer je manja disperzija

##############
# KORELACIJA #
##############
# U modelima koje cemo kasnije izucavati bitan je pojam korelacije, jer se u tim modelima pretpostavlja da sadasnja vrednost zavisi od vrednosti iz proslosti
# Ako otkrijemo tu korelaciju, to nam moze pomoci u prognoziranju i boljem uklapanju modela
# Ponovicemo korelaciju dve slucajne velicine
# Kovarijacija dve slucajne velicine - mera njihove linearne zavisnosti
# formula za kovarijaciju i koeficijent korelacije dve sl vel, uzoracke ocene
Herald <- read.delim("D:/.../Herald.txt")
# koncentracije benzoapirena i ugljen monoksida u 16 uzoraka vazduha na Herald Square-u na Menhetnu, oba se nalaze u izduvnim gasovima automobila
attach(Herald)
x <- CO
y <- Benzoa
n <- length(x)
#rucno izracunata kovarijacija, ili pomocu funkcije, ili asimptotski nepristrasna kovarijacija
sum((x - mean(x))*(y - mean(y))) / (n - 1)
cov(x,y)
mean((x - mean(x))*(y - mean(y)))
plot(x,y)
# Koeficijent korelacije - deljenje proizvodom standarnih devijacija
cov(x,y)/(sd(x)*sd(y))
cor(x,y)
##########
# Definicija strogo i slabo stacionarnog procesa
# Ova pretpostavka je veoma jaka i tesko bi bilo ispitati da li je neka serija realizacija strogo stacionarnog sl procesa
# Ispituje se zato slaba stacionarnost
# Da li je kovarijacija (odnosno korelacije) r_t i r_{t-l} invarijantna u vremenu (po t)
# Da li zavisi samo od l?
# Poznato nam je da slabo stacionarni procesi imaju konstantnu funkciju srednje vrednosti
# da je korelaciona funkcija, funkcija jednog argumenta (lag na engleskom) koji je prirodan broj
# smatramo da je interval uzorkovanja jedna jedinica vremena koja moze i treba da bude nesto malo, npr 0.1 sec
# Odsad stacionarnim procesom zovemo slabo stacionarne procese

# formula za (auto)kovarijacionu funkciju stacionarnog procesa
# formula za (auto)korelacionu funkciju stacionarnog procesa
# formula za ocenu ACVF
# formula za ocenu ACF

# Kazemo da serija nije korelisana ako je autokorelaciona funkcija jednaka 0, svuda, osim u nuli tu je po def ona jednaka 1

wave <- read.csv("D:/.../wave.txt", sep="")
# U date su razlike u vrednosti nivoa vode u sredini bunara koje su posledica simulatora talasa
# merenja su vrsena na svakih 0.1s i poslednji trenutak je 39.7s o
attach(wave)
plot(ts(waveht)) 
plot(ts(waveht[1:60]))
#sa grafika vidimo da nema fiksne frekvencije, ali da su uzastopne vrednosti relativno slicne
acf(waveht)
# dobijamo grafik autokorelacione funkcije - korelogram
# na x osi su prirodni brojevi
# na y osi su vrednosti autokorelacione funkcije sa datom zadrskom
# mozemo dobiti i odgovarajuci niz brojeva:
acf(waveht)$acf
# vrednost autokorelacione funkcije u 0 je 1
acf(waveht)$acf[1]
# vrednost autokovarijacione funkcije dobijemo kada funkciji acf dodamo jos jedan argument type
acf(waveht, type = "covariance")$acf[2]
# kao sto smo rekli, vrednost koeficijenta korelacije je mera linearne zavisnosti
# osim vrednosti same funkcije, mozemo nacrtati grafike da bismo se uverili u zavisnost
# linearnu, ili mozda neku drugu, izmedju vrednosti na rastojanju k vremenskih trenutaka
plot(waveht[1:396],waveht[2:397]) #ovo je sa zadrskom(lag) 1, na primer

## OSOBINE KORELOGRAMA ##
# 1 SAMO POLINOMIJALNI TREND
# linearni
y <- 1:100
y.ts<-ts(y)
acf(y.ts)
# polinom drugog stepena
y <- -y^2+y+5
y
acf(y)
# Vidimo da su sve korelacije pozitivne, statisticki znacajne i sporo opadaju
# Objasnjenje kao malopre jer je ovo opadajuca funkcija
# Generalno za strogo monotone tako izgleda korelogram

# 2 PERIODICNE-trigonometrijske
y <- seq(from=0,to=100,by=1)
acf(sin(y))
# oblik korelograma- kao cos

# 3 PERIODICNE- ponavljajuci niz brojeva
y <- rep(1:10,40)
acf(y)
# ponavlja se 10 brojeva, najvece su korelacije sa korakom 10*k

# Ispitajmo ovo na seriji AP u kojoj je bilo i trenda i sezonske varijacije
# Ona definitivno nije stacionarna, a videcemo kako trend i sezonska varijacija uticu na korelogram
data(AirPassengers)
AP <- AirPassengers
AP.ts <- ts(AP,freq=12,start=1958)
plot(AP.ts)
acf(AP.ts)
# sporo opadanje, zbog trenda
# izdvaja se autokorelacija sa korakom 12, zbog sezonske komponente
n <- (length(AP.ts)-6)
AP.decom <- decompose(AP.ts,type="mult")
acf(AP.ts[7:n]/AP.decom$trend[7:n],main="")
# Kada smo sklonili samo trend, a ostavili sezonsku komponentu, vidmo da je najveca pozitivna korelacija sa korakom 12
acf(AP.decom$random[7:n])
# i ovde vidimo kosinusni oblik na korelogramu sto moze da znaci da nam je los model
# Idealno je da serija greske bude beli sum, ali moze da znaci i da je model procesa greska AR(2) model (kasnije)
sd(AP[7:138])
sd(AP[7:138] - AP.decom$trend[7:138])
sd(AP.decom$random[7:138])
# izdvajanje deterministickih komponenti smanjuje disperziju pocetnog procesa
# posto je disperzija slucajne komponente mala, mozemo da zakljucimo da je model dobar
# i da je dobijeni oblik korelograma posledica AR(2) modela za proces greske

Fontdsdt <- read.csv("D:/.../Fontdsdt.txt", sep="")
# Posmatra se protok vode (m^3/s) (prosecno, mesecna serija, od 1.1909 do 12.1980)
# prethodno je ocenjen trend i sezonska komponenta, a adflow predstavlja samo slucajnu komponentnu
attach(Fontdsdt)
plot(ts(adflow), ylab = 'adflow')
acf(adflow, xlab = 'lag (months)', main="")
acf(w) # w - Gausov beli sum N(0,1)
# Jedino je statisticki znacajna korelacija sa korakom 1, sto se moze objasniti:
# u zavisnosti od podzemnih voda, cije se ponasanje (tj to koliko doprinose nasem jezeru) sporo menja
# pa ako su tokom uzastopnih meseci slicne vrednosti.
# Drugi nacin na koji se puni jezero su padavine, a padavine u uzastopnim mesecima nisu korelisane
# Ovakav oblik korelograma (eksponencijalno opadanje) je svojstveno za AR(1) model