#### Subjektive Priori für Präzision der Normalverteilung ####

#### Aufgabe b) ####

# Daten einlesen und plotten
wkurs <- read.table("Daten/r-jy-daily.dat")$V1
plot.ts(wkurs)
mean(wkurs)


# ML-Schätzer der ersten n0 Tage
n0 <- 10
kappaML <- n0/sum(wkurs[1:n0]^2)

# Priori-Parameter
a <- n0/2
b <- (a-1)/kappaML

# Priori zeichnen
s <- seq(0.1,3*1e5, 100)
plot(s, dgamma(s, shape = a, rate = b), type = "l",
     main = "Priori", xlab = expression(kappa), ylab = expression(p(kappa)))

# Posteriori-Parameter basierend auf übrigen Beobachtungen
wNew <- wkurs[-(1:n0)]
n <- length(wNew)
aPost <- a + n/2 
bPost <- b + sum(wNew^2)/2 # da mu = 0 angenommen wird

# Posteriori zeichnen
plot(s, dgamma(s, shape = aPost, rate = bPost), type = "l",
     main = "Posteriori", xlab = expression(kappa), ylab = expression(paste("p(", kappa , "| x)")))

# Relevanter Bereich
s <- seq(15000, 20000, 100)
plot(s, dgamma(s, shape = aPost, rate = bPost), type = "l",
     main = "Posteriori", xlab = expression(kappa), ylab = expression(paste("p(", kappa , "| x)")))


#### Aufgabe c) ####

# Schätzer für kappa
aPost/bPost # Posteriori-Erwartungswert
(aPost - 1)/bPost # alternativ: Posteriori-Median