#---------------------------------#
#   Additive und Multiplikative   #
#   Zufallsvariablen: Beispiel    #
#---------------------------------# 


# Wir wollen den Erwartungswert der s_n's 
# und der pi_n's mit Monte Carlo berechnen: 

n = 5                  # setzen wir dann noch auf 
                       # 10, 20, 50 oder 100..
N = 1000000            # 1 Mio MC-Simus, lassen wir fest 

phi = rnorm(n*N)
Phi = matrix(phi,N,n)
Phi2 = Phi^2
sn  = rowSums(Phi2)/n 
pin = apply(Phi2,1,prod)

sn
pin
summary(sn)
summary(pin)

mcsum_s  = cumsum(sn)/(1:N)
mcsum_pi = cumsum(pin)/(1:N)

Nplot = seq(from=100,to=N,by=100)     # 1 Mio Punkte sind etwas viel zum Plotten, 
length(Nplot)                         # wir nehmen nur jeden 100sten, sind dann 
                                      # nur 10000 Punkte
plot(mcsum_s[Nplot])
plot(mcsum_pi[Nplot])       # fuer n=5 klappt das noch..

# additional info:
plot(sort(sn)[Nplot])
plot(sort(pin)[Nplot])
tail(sort(pin),20)          # die groessten pin's..



# Wir plotten pin (und sn) als Funktion von n fuer  
# eine gegebene Realisation von Zufallszahlen und 
# machen das dann fuer, sagen wir, 20 Realisationen: 

n = 500 
for(i in 1:20)              # 20 Mal nacheinander
{
  phi = rnorm(n)
  phi2 = phi^2
  sn = cumsum(phi2)/(1:n)
  pin = cumprod(phi2) 

  # wir machen ein Plotfenster mit einem
  # 2x2 Bild-Array, das 4. Bild bleibt leer..
  par(mfrow=c(2,2))
  plot(sn,ylim=c(0,2))
  plot(pin,ylim=c(0,2))
  plot(1:n,log(pin))
  points(1:n,-1.27*(1:n),col="red") 
  readline("press enter to see next realization..")
}