#27.11.2000

library(ts)
library(SfS)

# Repetition vom letzten Mal:
# 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) #sollte jetzt installiert sein
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")

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)
}

#--------------------------------------------------------------------------
#Gleitende Mittel: y(t) = \sum coeff(u)*x(t-u). Wir plotten das Spektrum
#von y falls x weisses Rauschen ist, und vergleichen (x(t)) mit (y(t)),
#wenn x weisses Rauschen, bzw. eine harmonische Schwingung ist

f.transfer <- function(coef,nu)
{
  ## Purpose: Berechnung von \sum coeff(u)*exp(-i*2*pi*nu) 
  ##          (sogenannte Transferfunktion)
  ## -------------------------------------------------------------------------
  ## Arguments: coef=Koeffizienten, nu=Frequenz
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  2 Dec 96, 11:11
  p <- length(coef)-1
  z <- (cos(2*pi*nu)+sin(2*pi*nu)*1i)
  a  <- outer(z,0:p,FUN="^")
  a%*%coef}

nu <- seq(-0.5,0.5,0.01)
coef <- cos(0.3*pi*(0:7)) #Festlegung der Koeffizienten
coef <- coef/sqrt(sum(coef^2)) #Normierung der Koeffiziente
tr <- f.transfer(coef,nu)
plot(nu,Mod(tr)^2,type="l")
x <- ts(rnorm(500),start=1,freq=1)
y <- filter(x,filter=coef,sides=1)
plot(ts.union(x,y))
par(mfrow=c(1,2))
acf(x)
acf(y[8:500])
par(mfrow=c(1,1))
x <- ts(cos(2*pi*0.02*(1:100)),start=1,freq=1) #Modifizieren Sie 0.02 !

#Aliasing-effekt
time <- seq(0,8,0.01)
lambda <- runif(1)*pi
plot(time,cos(lambda*time),type="l")
points(0:8,cos(lambda*(0:8)))
lines(time,cos((lambda+2*pi)*time),lty=2)
lines(time,cos((lambda+4*pi)*time),lty=3)