f.kalman <- function(yt, model, f0, sm = F)
{
## Purpose: eindimensionaler Kalman-Filter
##          fuer AR(1)-Modell mit weissem Rauschen
## ----------------------------------------------------------------
## Arguments: yt: Beobachtungsreihe
##            model: Liste der Modellparameter
##            f0: Startwerte
##            sm: logische Variable, ob Glättung gerechnet wird
## ----------------------------------------------------------------
## Author: Markus Huerzeler, Date: Dec 96
  if(missing(model))
    stop("Modellspezifikation fehlt")
  if(missing(f0))
    stop("Startwerte fehlen")
  nn <- length(yt)
  ww <- numeric(nn)
  xfil <- c(f0$x0, ww)
  Rfil <- c(f0$R0, ww)
  Rpred <- ww
  xpred <- ww
  like <- 0
  for(i in 1:nn) {
    xpred[i] <- model$alpha * xfil[i]
    Rpred[i] <- model$sigma2 + (model$alpha)^2 * Rfil[i]
    xfil[i+1] <- xpred[i] + Rpred[i]/(Rpred[i] + model$tau2)*(yt[i]-xpred[i])
    ## ypred ist gleich xpred
    Rfil[i + 1] <- Rpred[i]/(Rpred[i] + model$tau2) * model$tau2
    like <- like - log(dnorm(yt[i], xpred[i], sqrt(Rpred[i] + model$tau2)))
  }
  xfil <- xfil[-1]
  Rfil <- Rfil[-1]
  res <- list(xfil=xfil, Rfil=Rfil, xpred=xpred, Rpred=Rpred, like=like)
  if(sm) {
    xsmooth <- xfil
    Rsmooth <- Rfil
    for(i in (nn - 1):1) {
      Jh <- (Rfil[i] * model$alpha)/Rpred[i+1]
      xsmooth[i] <- xfil[i] + Jh * (xsmooth[i+1] - xpred[i+1])
      Rsmooth[i] <- Rfil[i] + Jh * Jh * (Rsmooth[i+1] - Rpred[i+1])
    }
    res <- c(res, list(xsmooth=xsmooth, Rsmooth=Rsmooth))
  }
  return(res)
}

p.ts <- function(x, nrplots = max(1, min(8, n%/%400)), overlap = nk %/% 16,
                 main.tit = NULL, ylim = NULL, ylab = "",
                 quiet = FALSE, mgp = c(1.25, .5, 0), ...)
{
  ## Purpose: plot.ts with multi-plots + Auto-Title -- currently all on 1 page
  ## -------------------------------------------------------------------------
  ## Arguments: x      : timeseries [ts,rts,its,cts] or numeric vector
  ##            nrplots: number of sub-plots    [DEFAULT: in {1..8}, ~= n/400]
  ##            overlap: how much should subsequent plots overlap [DEFAULT:...]
  ##
  ## Depends on   mult.fig()
  ##
  ## ---> help page  ?p.ts
  ##
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date:  1 Jul 1994; 18.Dec,1998.

  if(is.null(main.tit)) main.tit <- paste(deparse(substitute(x)))
  isMat <- is.matrix(x)
  n <- if(isMat) nrow(x) else length(x)
  if(nrplots == 1)
    plot.ts(x, ..., main = main.tit, ylim = ylim, ylab = ylab)
  else if(nrplots <= 0)
    return(nrplots)
  else { # nrplots >= 2 :
    if(n <= 1) stop("`x' must have at least two points!")
    if(!is.ts(x)) x <- as.ts(x)
    ##-  do.dates <- !is.null(class(x)) && class(x) == "cts"
    ##-  if(do.dates) x <- as.rts(x)# dates() as below fails [S+ 3.4]
    scal <- (end(x) - (t1 <- start(x)))/(n-1)
    nk <- n %/% nrplots
    if(is.null(ylim))
       ylim <- range(pretty(range(x, na.rm = TRUE)))
    ##    --------
    pp <- mult.fig(mfrow=c(nrplots,1), main = main.tit, quiet= TRUE,
                   mgp = mgp, marP = c(-1,-1,-2,0))
    on.exit(par(pp $ old.par))
    for(i in 1:nrplots) {
      i0  <- max(0, (-overlap + (i-1)*nk)-1)
      in1 <- min(n, i*nk + overlap)-1
      st <- t1 + scal*i0 ##; if(do.dates) st <- dates(st)
      en <- t1 + scal*in1##; if(do.dates) en <- dates(en)
      if(!quiet) cat(i," -- start= ", format(st), "; end  =", format(en),"\n")
      plot(window(x, start= st, end = en), ylim = ylim, ylab = ylab, ...,
           plot.type = "single") # plot.type : for plot.mts only
    }
  }
}

mult.fig <- function(nr.plots, mfrow, mfcol,
         marP = rep(0, 4), mgp = c(1.5, 0.6, 0),
         mar = marP + 0.1 + c(4,4,2,1), oma = c(0,0, tit.wid, 0),
         main = NULL, tit.wid = if (is.null(main)) 0 else 1 + 1.5*cex.main,
         quiet = .Device == "postscript",
         cex.main = par("cex.main"), line.main = cex.main - 1/2,
         col.main = par("col.main"),
         font.main = par("font.main"),
         ...)
{
  ## Purpose: 'MULTiple FIGures' incl. title and other good defaults
  ## -------------------------------------------------------------------------
  ## Arguments: -- Either ONE of the first 3 arguments --
  ###  =========> help(mult.fig)
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, 1990 (UW, Seattle) -- 1995
  ## -------------------------------------------------------------------------

  use.row <- missing(mfcol)
  if (use.row)
    if (missing(mfrow)) {
      if (missing(nr.plots))
        stop("must either specify 'nr.plots', 'mfrow' or 'mfcol' !")
      else  mfrow <- n2mfrow (nr.plots)
    }
  old.par <<-
    if(use.row) par(mfrow = mfrow, oma= oma, mar = mar, mgp= mgp)
    else        par(mfcol = mfcol, oma= oma, mar = mar, mgp= mgp)
  if(!quiet) cat("Execute\n\t par(old.par) \n later to restore graphical par\n")
  ##---- now go ahead :
  if(!is.R())
      frame()
  if (!is.null(main)) {# Do title *before* first plot!
      if(is.R()) plot.new()
      mtext(main, side = 3, outer = TRUE,
            line = line.main,
            cex = cex.main,
            font = font.main, col = col.main, ...)
      if(is.R()) par(new=TRUE)# reverse `plot.new()' above
  }
  invisible(list(new.par = par(c("mfrow","mfcol","oma","mar","mgp")),
                 old.par = old.par))
}

f.a.saisonal <- function(nu)
{
  ## Purpose: Transferfunktion fuer Jahresfilter
  ## -------------------------------------------------------------------------
  ## Arguments: nu = frequenz
  ## -------------------------------------------------------------------------
  ## Author: Roberto Invernizzi, Date: 15 Dec 1998, 14:43
  j <- 1:5 ; h <- outer(j,2*pi*nu, "*"); mat <- matrix(cos(h),nrow=5)
  cos.sum <- apply(mat, 2, sum)
  1/12 + 1/12*cos(12*pi*nu) + 1/6*cos.sum
  }