## Programs for Applied Time Series Analysis Course

source('data.r')

fig.flood <- function(){
  pdf('flood-data.pdf', 
      width=15/2.54, height=15/2.54,
      family='GB1')
  on.exit(dev.off())
  opar <- par(mar=c(3,3,2,1), mgp=c(1.5,0.5,0))
  plot(flood.area1, lty=1,
       ylim=c(0,400),
       axes=FALSE,
       main='洪水灾害数据',
       xlab='年', ylab='面积(万亩)')
  box()
  axis(2)
  axis(1, at=seq(1948,1964,by=2))
  lines(flood.area2, lty=2)
  legend(1959, 380, lty=c(1,2), legend=c('受灾面积', '成灾面积'))
}

demo.flood <- function(){
  opar <- par(mar=c(3,3,3,1), mgp=c(1.5,0.5,0))
  on.exit(par(opar))
  plot(flood.area1, lty=1,
       ylim=c(0,400),
       axes=FALSE,
       main='洪水灾害数据',
       xlab='年', ylab='面积(万亩)')
  box()
  axis(2)
  axis(1, at=seq(1948,1964,by=2))
  lines(flood.area2, lty=2)
  legend(1959, 380, lty=c(1,2), legend=c('受灾面积', '成灾面积'))
}

## recode a factor using 0-1 matrix
factor.matrix <- function(y){
  y <- factor(y)
  n <- length(y)
  K <- length(levels(y))
  r <- matrix(0, nrow=n, ncol=K)
  r[cbind(seq(n), as.numeric(y))] <- 1
  r
}

demo.coal.data <- function(make.pdf=FALSE){
  if(make.pdf){
    pdf('coal-data.pdf', 
        width=15/2.54, height=15/2.54,
        family='GB1')
    opar <- par(mar=c(3,3,3,1), mgp=c(1.5,0.5,0))
    on.exit(dev.off())
  } else {
    opar <- par(mar=c(3,3,3,1), mgp=c(1.5,0.5,0))
    on.exit(par(opar))
  }
  y <- coal.consumption
  plot(y, lty=1, type='b',
       main='居民季度用煤消耗量',
       xlab='年', ylab='消耗量(吨)')
}

# 线性趋势加固定点周期项回归的程序
decomp.det1 <- function(y, log.trans=F){
  y0 <- y
  if(log.trans){
    y <- log(y0)
  }
  
  n <- length(y)
  s <- frequency(y)
  x1 <- 1:n
  if(s>1){
    x2 <- factor.matrix(cycle(y))
    x <- cbind(x1, x2[,-s])
    colnames(x) <- c("t", paste("s", 1:(s-1), sep=""))
  } else {
    x <- cbind(x1)
    colnames(x) <- "t"
  }
  
  lm.res <- lm(c(y) ~ x)
  y.fitted <- predict(lm.res)
  tr <- coef(lm.res)[1] + coef(lm.res)[2]*x1
  if(log.trans){
    y.fitted <- exp(y.fitted)
    tr <- exp(tr)
  }
  y.fitted <- ts(y.fitted, start=start(y), frequency=frequency(y))
  tr <- ts(tr, start=start(y), frequency=frequency(y))
  
  print(lm.res)
  
  plot(y0, ylab="")
  lines(y.fitted, col="blue")
  lines(tr, col="red")
  
  list(y=y0, y.fitted=y.fitted, coefs=lm.res)
}

demo.coal <- function(){
  oldpar <- par(mfrow=c(2,1),
                mar=c(3,3,3,1),
                mgp=c(1.5, 0.5, 0),
                oma=c(0,0,2,0))
  on.exit(par(oldpar))
  y <- coal.consumption
  ymore <- ts(c(y, rep(NA,4)), start=start(y), frequency=4)
  ymat <- matrix(c(y), byrow=TRUE, nrow=6, ncol=4)

  cols <- rainbow(20)
  ic <- 1

  ## 用同季度的值平均得到四个季节项
  get.season <- function(yd){ # input: Detrended series
    ymat <- matrix(c(yd), byrow=TRUE, ncol=4)

    ## season
    seas0 <- apply(ymat, 2, mean, na.rm=TRUE)

    seas0
  }

  ## 画去除趋势后的序列、季节项和随机项
  plot.season <- function(yd, seas0){ # input: Detrended series
    ## season
    seas <- rep(seas0, 6)
    seas <- ts(seas, start=c(1991,1), frequency=4)

    ## Error
    r <- c(yd) - seas
    r <- ts(r, frequency=4, start=c(1991,1))

    plot(yd, type='b', lwd=2,
         main="Detrended Series(black), Seasonal(red) and Random(cyan)",
         xlim=c(1991,1998))
    abline(v=1991:1996, col="gray")
    abline(h=0, col="gray")
    lines(seas, type="l", col="red")
    lines(r, type="l", col="cyan")
  }

  
  ## 用每年的平均值作为趋势项估计, 季度平均作为季节项
  tr0 <- apply(ymat, 1, mean, na.rm=TRUE)
  tr <- rep(tr0, each=4)
  tr <- ts(tr, frequency=4, start=c(1991,1))
  y.detrended <- y - tr

  seas0 <- get.season(y.detrended)

  tr.more <- ts(c(tr, rep(tr[length(tr)],4)),
                start=start(y), frequency=4)
  seas.more <- ts(rep(seas0, 7),
                  start=start(y), frequency=4)
  y.pred <- tr.more + seas.more

  plot(ymore, main="Series(black), Year average trend(red), Fit(green)",
       lwd=2,
       type="b", col="black")
  abline(v=1991:1997, col="gray")
  lines(tr.more, col="red", type="l")
  lines(y.pred, col="green", type="l")

  plot.season(y.detrended, seas0)
  mtext("每年平均值趋势，季度平均季节项", outer=TRUE)
  locator(1)


  ## 用时间的线性回归作为趋势项估计, 季度平均作为季节项
  ## linear trend method.
  yy <- c(y)
  tt <- seq(length(y))
  lmr <- lm(yy ~ tt)
  tr.more <- ts(predict(lmr, new.data=list(tt=seq(length(y)+4))),
                frequency=4, start=c(1991,1))

  ## season
  y.detrended <- y - tr.more[1:length(y)]
  seas0 <- get.season(y.detrended)
  seas.more <- ts(rep(seas0, 7),
                  start=start(y), frequency=4)
  y.pred <- tr.more + seas.more

  plot(ymore, main="Series(black), Linear trend(red) and Fit(green)",
       lwd=2,
       type="b", col="black")
  lines(tr.more, col="red", type="l")
  lines(y.pred, col="green", type="l")

  plot.season(y.detrended, seas0)
  mtext("线性回归趋势，季度平均季节项", outer=TRUE)
  locator(1)

  ## 二次曲线趋势
  yy <- c(y)
  tt <- seq(length(y))
  lmr <- lm(yy ~ tt + I(tt^2))
  tr.more <- ts(predict(lmr, new.data=list(tt=seq(length(y)+4))),
                frequency=4, start=c(1991,1))

  ## season
  y.detrended <- y - tr.more[1:length(y)]
  seas0 <- get.season(y.detrended)
  seas.more <- ts(rep(seas0, 7),
                  start=start(y), frequency=4)
  y.pred <- tr.more + seas.more

  plot(ymore, main="Series(black), Quadratic trend(red) and Fit(cyan)",
       lwd=2,
       type="b", col="black")
  lines(tr.more, col="red", type="l")
  lines(y.pred, col="green", type="l")

  plot.season(y.detrended, seas0)
  mtext("二次曲线回归趋势，季度平均季节项", outer=TRUE)
  locator(1)

  ## 加权移动平均趋势, c(1/8, 1/4, 1/4, 1/4, 1/8) filter
  yy <- c(y)
  tr <- ts(filter(yy, c(1/8, 1/4, 1/4, 1/4, 1/8), method="convolution"),
           frequency=4, start=c(1991,1))
  ind2 <- max(which(!is.na(tr)))
  tr.more <- ts(c(tr, rep(NA,4)),
                start=start(y), frequency=4)
  tr.more[(ind2+1):length(tr.more)] <- tr.more[ind2]
  y.detrended <- y - tr.more[1:length(y)]
  seas0 <- get.season(y.detrended)
  seas.more <- ts(rep(seas0, 7),
                  start=start(y), frequency=4)
  y.pred <- tr.more + seas.more

  plot(ymore, main="Series(black), MA trend(red), Fit(green)",
       lwd=2,
       type="b", col="black")
  abline(v=1991:1997, col="gray")
  lines(tr.more, col="red", type="l")
  lines(y.pred, col="green", type="l")

  plot.season(y.detrended, seas0)
  mtext("加权移动平均趋势，季度平均季节项", outer=TRUE)
}

## Poisson WN.
demo.pois <- function(){
  lambda <- 2.0
  n <- 60
  y <- ts(rpois(n, lambda) - lambda)
  plot(y, main="Poisson White Noise")
  invisible(y)
}


## Normal WN.
demo.norm <- function(){
  n <- 60
  y <- ts(rnorm(n))
  plot(y, main="Gaussian White Noise")
  invisible(y)
}

## Random phase harmonic WN
demo.harmonic <- function(){
  n <- 60
  b <- sqrt(2.0)
  aa <- c(4.0, 1.0, 0.5, 0.1)
  tt <- seq(n)
  for(a in aa){
    y <- ts( b * cos(a*tt + 2*pi*runif(n)) )
    plot(y, main=paste("Random harmonic WN: a=",
              a, sep=""))
    locator(1)
  }
  invisible(y)
}

## rectangle window MA filtering for harmonic + noise
demo.mafilt <- function(b=3, M=3){
  n <- 100
  ##om <- pi/7
  om <- pi/12
  sigma <- 1.0
  eps <- rnorm(n, 0, sigma)
  sn0 <- b^2 / (2*sigma^2)
  tt <- seq(n)
  u <- runif(1)
  signal <- b * cos(om*tt + 2*pi*u)
  y <- signal + eps
  filt <- rep(1/(2*M+1), 2*M+1)
  yf <- filter(y, filt, method="convolution")
  rg <- range(c(y, yf))
  sn <- round(b^2 / (2*sigma^2), 3)
  plot(tt, y, main=paste("MA filter: SN=", sn, sep=""),
       type="l", xlab='time', ylab='y',
       ylim=c(-b*1.5,b*1.5))
  lines(tt, signal, col="green", lwd=2)
  lines(tt, yf, col="red", lwd=2)
  legend(90, b*1.5, lty=c(1,1,1), lwd=c(1,2,2),
         col=c("black", "green", "red"),
         legend=c("序列", "信号", "滤波"))
}

## 频域观点：低频与高频
demo.freq.domain <- function(){
  oldpar <- par(mfrow=c(2,2)); on.exit(par(oldpar))
  m <- 10
  n <- 100+m+1
  fil <- rep(1/(m+1), m+1)
  eps <- rnorm(n)
  x <- filter(eps, fil, method="convolution", sides=1)
  eps <- eps[-seq(m+1)]
  x <- x[-seq(m+1)]
  ts.plot(x, col="black", main="Low frequency t.s.",
          xlab="time", ylab="y")

  fr <- seq(0, pi, length=100)
  sp <- numeric(length(fr))
  for(j in seq(along=fr)){
    sp[j] <- abs(sum(complex(mod=1, arg=-(0:m)*fr[j])))^2
  }
  sp <- sp/(2*pi*(m+1)^2)
  plot(fr, sp, type='l', main="Spectral density of LOW frequency t.s.",
       xlab="Frequency", ylab="Density")

  ## high frequency
  r0 <- complex(mod=1.1, arg=4/5*pi)
  poly1 <- c(1, -2*Re(r0), Mod(r0)^2)
  poly2 <- poly1 / poly1[3]
  poly2 <- rev(poly2)
  print(poly2)
  eps <- rnorm(n+100)
  y <- numeric(length(eps))
  for(j in 3:length(y)){
    y[j] = -poly2[2]*y[j-1] - poly2[3]*y[j-2] + eps[j]
  }
  y <- y[-seq(100)]
  ts.plot(y, main="High frequency t.s.",
          xlab="time", ylab="y")
  sp2 <- numeric(length(fr))
  for(j in seq(along=fr)){
    sp2[j] <- abs(sum(poly2 * complex(mod=1, arg=(0:2)*fr[j])))^2
  }
  sp2 <- 1/(2*pi) / sp2
  plot(fr, sp2, type='l', main="Spectral density of HIGH frequency t.s.",
       xlab="Frequency", ylab="Density")
}

demo.flute <- function(){
  require(seewave)
  require(sound)
  ##opar <- par(mfrow=c(2,2), mar=c(2,2,1,0.2), mgp=c(2,1,0))
  ##on.exit(par(opar))
  par(mfrow=c(1,1))
  
  ## 用Audacity打开后播放，播放三个片段。
  ## 显示频谱。
  mus <- loadSample("music/01music-demo2.wav")
  print(mus)

  f <- rate(mus)
  mi1 <- cutw(mus, from=0.125, to=1.125,
              Sample=TRUE, output="Wave")
  mi1a <- cutw(mus, from=0.500, to=0.5500,
               Sample=TRUE)
  so <- cutw(mus, from=1.190, to=1.420,
             Sample=TRUE, output="Wave")
  mi2 <- cutw(mus, from=1.640, to=2.660,
              Sample=TRUE, output="Wave")
  ## cutw 取出声音片段。结果为一列的矩阵（列向量）。

  oscillo(mi1a, f=f,
          title='0.05秒的声音时间序列图') ## 显示0.1秒的序列图。
  locator(1)
  spec.pgram(ts(mi1a, freq=f), spans=c(5,7),
             xlim=c(0,10000), main="用spec.pgram进行谱估计",
             xlab="")
  locator(1)
  ## 二维谱线图，动态显示随时间变化的主要谱峰位置
  spectro(mus, f=f, flim=c(0,2), main="动态谱峰图")
  locator(1)
  ##saveSample(mi1, filename="music/91mi1.wav")
  ## 直接生成一段单音，C，523Hz.
  f <- 8000
  y <- sin(2*pi*523*seq(0,1, length=f*1))
  mip <- repw(y, f=f, times=1, Sample=TRUE)
  par(mfrow=c(1,1))
  oscillo(cutw(mip, f=f, from=0, to=0.05, Sample=TRUE),
          f=f,
          title='直接用sin函数生成523Hz单音')

  locator(1)
  spec.pgram(ts(mip, freq=f), spans=c(5,7),
             xlim=c(0,2000), main="频率523Hz单音(C)的谱估计",
             xlab="")
  abline(v=523, col="red")
  locator(1)
  spectro(mip, f=f, flim=c(0,1.0), main="523Hz单音(C)的动态谱峰图")
  ##saveSample(mip, filename="music/91mi-pure.wav")

  ## 播放生成的单音。音调与mi1相同但有明显的电子味道。
}


demo.song <- function(){
  require(seewave)
  ##require(sound)
  opar <- par(mfrow=c(2,2), mar=c(2,2,1,0.2), mgp=c(2,1,0))
  on.exit(par(opar))

  freq1 <- 262
  freq1b <- 523
  freq2 <- 1568
  freq2b <- 784
  
  f <- 8000 # 采样频率
  dur <- 1  # 长度描述

  ## 直接生成一段单音，低音C，262Hz.
  y1 <- sin(2*pi*freq1*seq(0,1, length=f*dur))
  sa1 <- repw(y1, f=f, times=1, Sample=TRUE)
  oscillo(cutw(sa1, f=f, from=0, to=0.02, Sample=TRUE), f=f,
          title='低音C单音')
  spec.pgram(ts(sa1, freq=f), spans=c(5,7),
             xlim=c(0,4000), main="", xlab="")
  abline(v=freq1, col="red")
  ##saveSample(sa1, filename="music/92C262-pure.wav")


  ## 直接生成一段单音，高音G，1568Hz.
  y2 <- sin(2*pi*freq2*seq(0,1, length=f*dur))
  sa2 <- repw(y2, f=f, times=1, Sample=TRUE)
  oscillo(cutw(sa2, f=f, from=0, to=0.02, Sample=TRUE), f=f,
          title='高音G单音')
  spec.pgram(ts(sa2, freq=f), spans=c(5,7),
             xlim=c(0,4000), main="", xlab="")
  abline(v=freq2, col="red")
  ##saveSample(sa2, filename="music/92G1568-pure.wav")
  locator(1)

  ## 生成复音。
  y1b <- 0.5*sin(2*pi*freq1b*seq(0,1, length=f*dur))
  sa1b <- repw((y1+y1b)/2, f=f, times=1, Sample=TRUE)
  oscillo(cutw(sa1b, f=f, from=0, to=0.02, Sample=TRUE), f=f,
          title='低音C和中音C的复音')
  spec.pgram(ts(sa1b, freq=f), spans=c(5,7),
             xlim=c(0,4000), main="", xlab="")
  abline(v=c(freq1,freq1b), col="red")
  ##saveSample(sa1b, filename="music/92C262-523-pure.wav", overwrite=TRUE)
  
  ## 生成复音。
  y2b <- 0.5*sin(2*pi*freq2b*seq(0,1, length=f*dur))
  sa2b <- repw((y2+y2b)/2, f=f, times=1, Sample=TRUE)
  oscillo(cutw(sa2b, f=f, from=0, to=0.02, Sample=TRUE), f=f,
          title='高音G和中音G的复音')
  spec.pgram(ts(sa2b, freq=f), spans=c(5,7),
             xlim=c(0,4000), main="", xlab="")
  abline(v=c(freq2,freq2b), col="red")
  ##saveSample(sa2b, filename="music/92G1568-784-pure.wav", overwrite=TRUE)
  
  invisible()
}

### 1.8节的单个频率的调和平稳序列的谱函数图形。
dspec <- function(){
  pdf(file='181.pdf', width=20/2.54, height=12/2.54)
  opar <- par(mfrow=c(1,1),
              mar=c(2, 0, 0, 0), mgp=c(0,0,0))
  on.exit({par(opar);dev.off()})
  
  plot(c(-pi-0.5, pi+0.5), c(-0.1, 1.25),
       type="n", axes=F, main="", xlab="", ylab="")
  arrows(-pi-0.5, 0, pi+0.5, 0, lwd=2)
  arrows(0, -0.1, 0, 1.15, lwd=2)
  text(0.1, -0.02, "0")
  lines(c(-pi, -pi), c(0, 0.02))
  lines(c(pi, pi), c(0, 0.02))
  text(-pi, -0.03, expression(-pi))
  text(pi, -0.03, expression(pi))
  lam0 <- 2
  lines(c(-lam0, -lam0), c(0, 0.02))
  lines(c(lam0, lam0), c(0, 0.02))
  text(-lam0, -0.04, expression(-lambda[0]))
  text(lam0, -0.04, expression(lambda[0]))
  lines(c(-0.07,0.07), c(1,1))
  text(-0.2, 1.02, expression(sigma^2))
  lines(c(-pi, -lam0), c(0,0), lwd=3)
  points(-lam0, 0, pch=1)
  lines(c(-lam0, lam0), c(0.5,0.5), lwd=3)
  points(lam0, 0.5, pch=1)
  lines(c(lam0, pi), c(1,1), lwd=3)
  text(0.1, 1.2, expression(F(lambda)))
}
##dspec()

### 美国泛美航空公司国际航班订票数（单位：千名）
demo.airline.data <- function(){
  opar <- par(mar=c(3,3,3,1), mgp=c(1.5,0.5,0))
  on.exit(par(opar))

  data(AirPassengers)
  y <- AirPassengers
  cat("国际航班订票数（千人）：\n")
  print(y)
  plot(y, main='国际航班订票数（千人）',
       ylab='')
  ## 开始点
  cat("序列开始时间：\n")
  print(start(y))
  cat("序列结束时间：\n")
  print(end(y))
  cat("序列采样频率：\n")
  print(frequency(y))

  browser()
}

### Maine 1996.1--2006.8失业率
demo.Maine.data <- function(){
  opar <- par(mar=c(3,3,3,1), mgp=c(1.5,0.5,0))
  on.exit(par(opar))
  y <- Maine
  plot(y, main='Maine失业率', ylab='')
}

### 用国际航班订票数据分解后的随机项演示谱峰
demo.month.spec <- function(){
  opar <- par(mfrow=c(1,1),
              mar=c(3,4,3,0.5),
              mgp=c(2,0.5,0))
  on.exit(par(opar))
  y <- AirPassengers
  res <- decompose(y)
  r <- res$random
  plot(r, main="航班数据的随机项")
  locator(1)
  r2 <- r[!is.na(r)]
  spec.ar(r2, main="谱密度估计")
  abline(v=1/12)
  text(0.12, 30, expression(1/12))
  ## 说明：谱密度估计的纵轴用了对数刻度，
  ## 频率范围从0到0.5, 对应角频率的0到pi。
  ## 角频率除以2\pi变成频率。
}
##demo.month.spec()


demo01 <- function(){
  demos <- list("洪水灾害数据图"=demo.flood,
                "国际航班订票数据"=demo.airline.data,
                "Maine失业率数据"=demo.Maine.data,
                "煤炭消耗量数据"=demo.coal.data,
                "煤炭消耗量时间序列分解"=demo.coal,
                "Poisson WN."=demo.pois,
                "Gaussian WN."=demo.norm,
                "Harmonic WN."=demo.harmonic,
                "矩形窗滑动平均滤波"=demo.mafilt,
                "笛子声音样本和频谱" = demo.flute,
                "人为正弦波声音和频谱" = demo.song,
                "低频时间序列和高频时间序列" = demo.freq.domain,
                "月度数据的谱密度" = demo.month.spec
                )
  ndemos <- length(demos)
  prompts <- paste("   ", seq(ndemos), ". ", names(demos), "\n", sep="")
  cat("Demos of Chapter 1:\n",
      prompts)
  sel <- scan(n=1, quiet=TRUE)
  demos[[sel]]()
}