#20.11.2000
library(ts)
library(SfS)

#Repetition vom letzten Mal
#Autokorrelationen: Versuchen Sie, die Zusammenhaenge zwischen dem
#Zeitreihenplot und der Autokorrelationsfunktion zu verstehen. Die Funktion
#f.lagplot(x) zeigt den Zusammenhang zwischen x(t) und x(t+h) fuer kleine
#h=1,2,..,12, was je nach Zeitreihe informativer ist als die Autokorrelation.

#1: Luchs
par(mfrow=c(2,1))
data("lynx")
acf(lynx)
acf(log(lynx))
par(mfrow=c(1,1))
lag.plot(log(lynx),lags=12,layout=c(3,4))
#2: Airlines
airline <- ts(scan.url("http://stat.ethz.ch/Teaching/Datasets/airline.dat"),
              start=(1949,1),freq=12)
acf(log(airline))
acf(diff(diff(log(airline)),lag=12))
airl.stl <- stl(log(airline),s.window=7,t.window=21)
acf(airl.stl$time.series[,3])
#3: Ozeanwellen: "The change in the level of ambient noise in the ocean"
# (4 Messungen pro Sekunde), aus Percival& Walden, Spectral Analysis for
# Physical Applications
ocwave <- ts(scan.url("http://stat.ethz.ch/Teaching/Datasets/ocwave.dat"),
          start=1,deltat=0.25)
plot(ocwave)
acf(ocwave)
acf(ocwave,lag.max=200)
#4: AR
x <- filter(ts(rnorm(200)),filter=0.5,method="recursive",init=rnorm(1))
acf(x)             
#5: Arch
              f.arch.sim <- function(n=500,init=200,a)
{
  ## Purpose: Simulating the Arch(1) process
  ## -------------------------------------------------------------------------
  ## Arguments: n=length of the series, init=length of the initial segment
  ## that is discarded, a= Parameter of the model
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  5 Feb 97, 16:49
  m <- n+init
  x <- rnorm(m)
  for (i in (2:m)) x[i] <- sqrt(1+a*x[i-1]^2)*x[i]
  ts(x[(init+1):m],start=1)
}
x <- f.arch.sim(a=0.8)
acf(x)
acf(x^2)
lag.plot(x,lags=4,layout=c(2,2))          
#6: Lineare Zunahme
acf(1:100,lag.max=50)
   #Warum wird eine geschaetzte Autokorrelation negativ, obwohl stets x(t+h)
   # groesser als x(t) ?           

#7: Deterministisches chaotisches System
f.logist <- function(n,start=runif(1))
{
  ## Purpose: Computing the iterated logistic map x(t+1)=4*x(t)*(1-x(t))
  ## -------------------------------------------------------------------------
  ## Arguments: start=starting value, length=number of iterates
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 11 Nov 96, 14:59
  x <- rep(0,n)
  x[1] <- start
  for (i in (2:n)) {x[i] <- 4*x[i-1]*(1-x[i-1])}
  x
}          
x.det <- ts(f.logist(200),start=1,deltat=1)
plot(x.det)
acf(x.det)
lag.plot(x.det,lags=12,layout=c(3,4))
--------------------------------------------------------------------------
# Simulationen zur Untersuchung der Variabilitaet des arithmetischen
# Mittels: Die simulierten Reihen werden so standardisiert, dass Var(X_t)=1,
# und dann werden die arithmetischen Mittel mit der Laenge der Reihe hoch
# H-1 multipliziert. Man soll nachpruefen, ob die Streuungen der so
# standardisieren arithmetischen Mittel unabhängig von der Länge der Reihe
# wird, wie gross diese Streuungen sind, und ob genäherte Normalverteilung gilt.

f.sim.mean <- function(n=800,H=0.5,x.sim,...)
{
  ## Purpose: Simulating a time series of length n and computing the
  ## standardized mean for subseries of length n/8, n/4, n/2, n
  ## -------------------------------------------------------------------------
  ## Arguments: H=exponent for standardizing the mean
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 20 Nov 2000, 11:49
  x <- x.sim(n,...)
  x <- x/sqrt(var(x))
  k <- floor(n/8)
  m <- rep(0,4)
  for (i in (1:4)){
    m[i] <- mean(x[1:k])*k^(1-H)
    k <- 2*k}
  m}

#Simulation mit i.i.d. t-verteilten Daten (5 Freiheitsgrade)
resv <- matrix(0,nrow=200,ncol=4)
for (i in (1:200)) resv[i,] <- f.sim.mean(x.sim=rt,df=5)
boxplot(data.frame(resv))
par(mfrow=c(2,2))
qqnorm(resv[,1],main="n=100")
qqnorm(resv[,2],main="n=200")
qqnorm(resv[,3],main="n=400")
qqnorm(resv[,4],main="n=800")
par(mfrow=c(1,1))

#Simulation mit einem autoregressiven Prozess
f.ar.sim <- function(n, beta, sigma = 1.0)
{
    x <- ts(rnorm(n+100, 0, sigma), start = -99)
    x <- filter(x, beta, method = "recursive")
    as.ts(x[-(1:100)])
}

resv <- matrix(0,nrow=200,ncol=4)
for (i in (1:200)) resv[i,] <- f.sim.mean(x.sim=f.ar.sim,beta=0.8)
#Anschauen der Resultate wie vorher

#Simulation mit einem Prozess, bei dem R(k) asymptotisch abfaellt wie
#k^(2H-2). Die Funktion zur Simulation eines solchen Prozesses verwendet
#den Parameter d=H-0.5
library(fracdiff)
f.lrd.sim <- function(n,d) fracdiff.sim(n,ar=0,ma=0,d)$series


resv <- matrix(0,nrow=400,ncol=4)
for (i in (1:200)) resv[i,] <- f.sim.mean(x.sim=f.lrd.sim,d=0.3)
for (i in (1:200)) resv[i,] <- f.sim.mean(H=0.8,x.sim=f.lrd.sim,d=0.3)

#Simulationen zur Untersuchung des Verhaltens der geschaetzten
#Autokorrelationen 

f.sim.acf <- function(n=200,lag=10,x.sim,...)
{
  ## Purpose: Simulation of the autocorrelation function for various models
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 20 Nov 2000, 14:30
  x <- x.sim(n,...)
  acf(x,lag.max=lag,plot=F)$acf[-1]}

resv <- matrix(0,nrow=200,ncol=10)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=f.ar.sim,beta=0.8)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=f.arch.sim,a=0.8)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=rnorm)
boxplot(data.frame(resv))
plot(resv[,1],resv[,2],xlab="ac for lag 1",ylab="ac for lag 2")
qqnorm(resv[,1],main="lag=1")