###################
### CETVRTI CAS ###
###################
# Za ispitivanje osobina nekog modela korisno je koristiti se simulacijama
# Do sada smo imali realne serije, a u narednom kodu mi simuliramo serije i onda 
# ispitujemo njihove osobine (ovo je prakticni pristup, a ne teorijsko ispitivanje) i pokusavamo
# da ih uklopimo u model (iz kog smo prvobitno simulirali)

# Beli sum, Gausov beli sum
set.seed(1) # cemu sluzi?
w <- rnorm(100)
plot(w, type = "l")
# ovim se samo ilustruje varijabilnost uzorka - iako generisemo tacno iz N(0,1) raspodele
# histogram bas i ne odgovara gustini
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) #ili moze ovako da se doda gustina na histogram

# Kako prepoznati beli sum?
# funkcija srednje vrednosti/autokovarijaciona/autokorelaciona funkcija za beli sum
# pri simulacijama ove vrednosti nece biti kao teorijske, dovoljno je da budu dovoljno blizu
set.seed(2)
acf(rnorm(100))$acf
# iako smo simulirali tacno iz N(0,1) raspodele, kod korelograma dobijamo 
# 1 statisticki znacajnu korelaciju (k=7)
# po ovome trebalo bi odbaciti nultu hipotezu o tome da su sve korelacije 0
# napravila bi se greska prve vrste, cija je verovatnoca 0.05

# ako imamo seriju reziduala i verujemo da je beli sum, treba oceniti \sigma^2 
x <- rnorm(100,mean=0,sd=5)
var(x)

## Slucajno lutanje
x <- w <- rnorm(1000)
for (t in 2:1000) x[t] <- x[t-1] + w[t]
plot(x, type = "l") # rastuci trend
acf(x)$acf[2] # crta grafik i ispisuje vrednost r(1) - koja je bliska 1
# kod slucajnog lutanja korelacije polako opadaju od r(0)=1 
acf(x, type = "covariance")$acf[2]
# crtanje autokovarijacione funkcije i ispisivanje kovarijacije c(1)

# sledeci deo koda - implementirane funkcije acf i acvf
sample.acvf <- function(x,k){
  n <- length(x)
  s <- 0
  for(i in 1:(n-k)){
    s<-s+(x[i]-mean(x))*(x[i+k]-mean(x))
  }
  s<-s/n
  return(s)
}
sample.acvf(x,1)

sample.acf<-function(x,k){
  return(sample.acvf(x,k)/sample.acvf(x,0))
}
sample.acf(x,1)
# acf za slucajno lutanje ima odredjeni izgled - sporo opadajuce autokorelacije
# ali da li je neki proces random walk mozemo proveriti i koriscenjem diferencirane serije
# diferencirana serija dobija se sa: \Delta x_t = x_t-x_{t-1}
# Kod slucajnog lutanja diferencirana serija je bas beli sum, a znamo kako za njega izgleda korelogram
acf(diff(x))

## Braunovo kretanje
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)*sqrt(0.001)
}
plot(t,B, type="l")

# Primer 
# exchange rate - koliko jedne valute se dobija za jedinicnu vrednost druge
Z <- read.csv("D:/.../1cas/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")
# ukoliko prvo pretpostavimo da je odgovarajuci model slucajno lutanje, za proveru ispitujemo diferenciranu serija
acf(diff(Z.ts))
# Da li je diferencirana serija dobro modelovana ovim procesom, slicnim belom sumu?
# Potrebno je proveriti:
    # 1. Da li su nekorelisane vrednosti procesa?
    # 2. Da li su normalno raspodeljene?
    # 3. Da li je disperzija konstantna tokom vremena (jednaka \sigma^2)?
# Kako proveravamo ove osobine?
    # 1. na osnovu korelograma mozemo zakljuciti da li ili ne postoji serijska korelacija
    # 2. najjednostavnije je graficko testiranje - qqnorm grafika
    # 3. ovim se necemo baviti dublje, ako se na osnovu grafika reziduala moze primetiti da nema
    # razlike u variranju tokom vremena, zakljucicemo da je disperzija konstantna tokom vremena
# Ipak, trazicemo bolji model Holt-Winters metodom
Z.hw <- HoltWinters(Z.ts, gamma = FALSE) 
# ima smisla sto smo stavili gamma=FALSE jer se sa grafika ne uocava sezonska komponenta
Z.hw$alpha 
Z.hw$beta
# dobija se alpha=1
# ovo je kao BM samo sa dodatnim clanom
# raspisati jednacine modela za proveru
acf(resid(Z.hw))
# sada se dobija i bolji korelogram - nijedna korelacija nije statisticki znacajna

## Slucajno lutanje sa driftom
#Primer 2
HP <- read.csv("D:/.../HP.txt", sep="")
attach(HP)
# cene akcija za Hewlett-Packard Company (do 7. juna2007)
plot (as.ts(Price))
# Diferenciranjem ove serije dobijamo proces koji se od belog suma razlikuje samo po tome sto ima srednju vrednost \delta, a ne 0.
DP <- diff(Price)
# Provere za beli sum:
acf(DP)
plot(as.ts(DP))
mean(DP) # ocena drifta
qqnorm(DP,lty=2) #jedino ovo bas i ne prolazi
mean(DP) + c(-2, 2) * sd(DP)/sqrt(length(DP)) # 95% interval poverenja za drift delta
# pod pretpostavkom da DP ima normalnu raspodelu
# interval poverenja u nasem slucaju je skup tacaka x za koje ne moze da se odbije nulta hipoteza 
# da je srednja vrednost jednaka x, 
# dobija se interval koji ne sadrzi 0, sto znaci da se za tacku 0 hipoteza se odbija 
# dakle, odbacujemo hipotezu da je drift=0