#### Bayesianische Kredibilitätsintervalle für Präzision der Normalverteilung ####

#### Vorarbeiten ####

# 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

# 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 b) ####

# R-Skript mit Hilfsfunktion zur Berechnung des HPD-Intervalls laden
source("....R")

# 95-% HPD-Intervall für Posteriori
HPD <- gammaHPD(0.95, aPost, bPost)
HPD
lines(HPD$interval, c(0,0), lwd = 5, col = 4) # plotten

# Vergleich mit gleichendigem 95%-Intervall
u <- 0.05/2
Gleich <- c(qgamma(u, shape = aPost, rate = bPost), qgamma(0.95 + u, shape = aPost, rate = bPost))
Gleich # Grenzen
diff(Gleich) # Intervall-Länge
lines(Gleich, c(1e-5,1e-5), lwd = 5, col = 3) # plotten
legend("topright", c("HPD", "Gleichendig"), col = c(4,3), lwd = 5)