## The Stock Market Data (ISLR)

###
library(ISLR2)
?Smarket

head(Smarket)
## Direction=Up akko Today>0
dim(Smarket)
summary(Smarket)

###
cor(Smarket[, -9])
# velika korelacija izmedu Volume i Year; 
# mala korelacija Today i Lag_k (ocekivano)

boxplot(Volume~Year, data=Smarket)
### Volume = Volume of shares 
#   traded (number of daily shares traded in billions) 


attach(Smarket)
### podatke prije 2005 koristimo za trening
train <- (Year < 2005)
Smarket.2005 <- Smarket[!train, ]
dim(Smarket.2005)
# 252 povrata u 2005
n=1250-252 ## velicina skupa za trening
Direction.2005 <- Direction[!train]

############################ logisticka regresija ###################################

### Neka je Y=Direction te pretpostavljamo model 
##  P(Y="Up" | Lag1,...., Lag5, Volume) 
##  = h(b0 + b1*Lag1 + ... + b5*Lag5+b6*Volume), za
##  h(z) = exp(z)/(1+exp(z)) (inverz logit funkcije)

glm.fits <- glm(
  Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5+Volume,
  data = Smarket, family = binomial, subset = train
)


summary(glm.fits)
### b_lag1<0 implicira da se vjerojatnost za "up" smanjuje
## ako je povrat u prethodnom danu bio (jako) pozitivan

### p-vrijednost u zadnjem stupcu odnosi se na 
##  testiranje H_0: bj=0 pri cemu je testna
##  statistika slicna kao kod lin. regresije 
##  (tzv. Waldova statistika), te rezultate slicno intepretiramo
##  -> niti jedna se zasebno ne cini znacajna 
##  (najmanja p-vrijednost za Lag1)

### Ipak, ne mozemo odmah zakljuciti da niti 
##  jedna kovarijata nije
##  stat. znacajna u modeliranju odziva u ovom modelu


## tzv. nul-model
glm.fits0 <- glm(
  Direction ~ 1,
  data = Smarket, family = binomial, subset = train
)

anova(glm.fits0,glm.fits, test="Chisq")
### smanjenje devijance je jako malo (2.1601),
##  tj. p-vrijednost jako velika, pa
##  ne mozemo odbaciti nultu hipotezu da je 
##  model samo sa slobodnim clanom
## (dakle, bez kovarijata) dovoljan 

### pogledajmo ipak kakva je prediktivna sposobnost modela

## procjena greske na testnom skupu
glm.probs <- predict(glm.fits, Smarket.2005,
                     type = "response")

contrasts(Direction)
### dakle 1="up"

glm.pred <- rep("Down", 252)
### ako koristimo 0-1 gubitak, onda
##  za predikcije gledamo je li 
##  P("up" mid X=x) veća ili manja od 0.5
glm.pred[glm.probs > .5] <- "Up"

table(glm.pred, Direction.2005) 
### kontingencijska tablica za vektor (glm.pred, Direction.2005)
#          Direction.2005
# glm.pred Down  Up
#   Down   77   97
#   Up     34   44
### procjena testne greske P(Y_hat != Y) 
## (dakle, uz 0-1 funkciju gubitka) je
mean(glm.pred != Direction.2005)
##  0.5198413 = (34+97)/252 --> losije od slucajnog pogadanja!


### probajmo:
glm.fits <- glm(Direction ~ Lag1 + Lag2, data = Smarket,
                family = binomial, subset = train)
glm.probs <- predict(glm.fits, Smarket.2005,
                     type = "response")
glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"
table(glm.pred, Direction.2005)
#          Direction.2005
# glm.pred Down  Up
# Down       35  35
# Up         76 106

mean(glm.pred != Direction.2005)
### [1] 0.4404762 --> izgleda puno bolje, ali
## kada bi stavili f_hat(x):= "up" za sve x, 
## dobili bi istz testnu gresku (to je vjerojatno slucajnost)
mean(Direction.2005 != "Up")
# [1] 0.4404762

## Ipak, uocimo da je u nasem modelu procjena 
## za uvjetnu gresku P(Y != Y_hat | Y_hat="Up")
76/(76+106)
# [1] 0.4175824

## Primjer predikcije

h=function(z) exp(z)/(1+exp(z))

# nas model je sada 
# P(Y="Up" | Lag1, Lag2) = h(b0 + b1*Lag1 + b2*Lag2)

b=coef(glm.fits)
# procjena za P(Y="Up" | Lag1=1.2, Lag2=1.1) je
h(b[1] + b[2]*1.2 + b[3]*1.1)
# 0.4791462 

# ili preko predict:
predict(glm.fits,
        newdata =
          data.frame(Lag1 = 1.2,  Lag2 = 1.1),
        type = "response"
)
# 0.4791462 


######################################### LDA ###############################

library(MASS)
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket,
               subset = train)

lda.fit
# dakle, ako je 0="Down", 1="Up", 
# 1. pi0_hat=0.491984, pi1_hat=0.508016
# 2. mi0_hat=(0.04279022, 0.03389409) approx E(X | Y=Down), 
#    mi1_hat=(-0.03954635, -0.03132544) approx E(X | Y=Up)
### ima li 2 smisla?

# plot(lda.fit)

lda.pred <- predict(lda.fit, Smarket.2005)
names(lda.pred)

# lda.pred$posterior daje 
# P(Y= "Down" | X=x) i P(Y= "Up" | X=x)

lda.class <- lda.pred$class
table(lda.class, Direction.2005)

mean(lda.class != Direction.2005)
# [1] 0.4404762 -- isto kao kod log. regresije jer 
# daju potpuno iste predikcije na testnom skupu:

table(lda.class, glm.pred)


########################## QDA ########################################
qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket,
               subset = train)
qda.fit

qda.class <- predict(qda.fit, Smarket.2005)$class
table(qda.class, Direction.2005)

mean(qda.class != Direction.2005)
# [1] 0.4007937 -- bolje nego log. regr. i LDA

######################## Naive Bayes #####################
library(e1071)
nb.fit <- naiveBayes(Direction ~ Lag1 + Lag2, data = Smarket,
                     subset = train)
### pretpostavlja se da su uvjetno na Y=k, 
##  kovarijate nezavisne i normalno distribuirane
nb.fit
### npr. uvjetno na Y="Down", 
## procijenjeno ocekivanje i varijanca za 
## Lag1 su 0.04279022 i 1.227446.

nb.class <- predict(nb.fit, Smarket.2005)
table(nb.class, Direction.2005)

mean(nb.class != Direction.2005)
## [1] 0.4087302 -- samo malo losije od QDA, uocimo:
table(qda.class, nb.class)

## izgleda da LDA i log. regresija u ovom slucaju 
## imaju preveliku pristranost jer granica odluke 
## nije priblizno linearna

## procjena vjerojatnosti P(Y=k | X=x) na testnom skupu
nb.preds <- predict(nb.fit, Smarket.2005, type = "raw")
nb.preds[1:5, ]


#######################  kNN ################################
### slicno kao kod regresije

#sintaksa je malo drugacija

library(class)
train.X <- cbind(Lag1, Lag2)[train, ]
test.X <- cbind(Lag1, Lag2)[!train, ]
train.Direction <- Direction[train]

set.seed(1) ## kako bi mogli reproducirati rezultate jer
# u slucaju da ima vise podatkaka koji su k-ti susjed, 
# slucajno se odabire jedan
knn.pred <- knn(train.X, test.X, train.Direction, k = 1)
table(knn.pred, Direction.2005)

mean(knn.pred != Direction.2005)
# [1] 0.5 -- k=1 daje preveliku varijancu

set.seed(5)
knn.pred=numeric(300)
for (i in 1:300) {
  pred <- knn(train.X, test.X, train.Direction, k = i)
  knn.pred[i]=mean(pred != Direction.2005)
}
c(which.min(knn.pred), min(knn.pred))

### [1] 242.0000000   0.3888889 -- izgleda super! 
## U cemu je problem?

### -> Ovo nije dobra procjena za testnu gresku!
### Opcenito, skup za odabir parametra kompleksnosti 
##  i procjenu greske ne smije biti isti! 
### Koristimo dodatni skup za tzv. validaciju
### DZ -- trebalo bi k odabrati unakrsnom validacijom
## na trening skupu, te onda procijeniti gresku na testnom skupu