Untitled

y <- c(AirPassengers) # 月間航空旅客数
n <- length(y)        # 月数
n.tr <- 120           # 訓練月数
n.te <- n - n.tr      # テスト月数
t <- 1:n              # 月
ii.tr <- 1:120        # 訓練期間インデックス
ii.te <- 121:n        # テスト期間インデックス
matplot(y, type = 'o', pch = 1, col = 2)

fq <- 12

library(MASS)
bc <- boxcox(y ~ t, data = data.frame(y = y[ii.tr], t = t[ii.tr]))

lambda <- round(bc$x[which.max(bc$y)], 1)
lambda

## [1] 0

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

y.bc <- BoxCox(y[ii.tr], lambda)

matplot(y.bc, type = 'o', pch = 1, col = 2, ylab = 'y')

d <- 1

dy <- diff(y.bc, diff = d)

matplot(dy, type = 'o', pch = 1, col = 2, ylab = 'y')

ys <- ts(10 * sin(2 * pi * t / 10) + rnorm(n))

dir.create('acf_lag', showWarnings = F)

for (k in 1:25)
{
  png(paste0('acf_lag/acf_lag', k, '.png'))

  ys.lag <- lag(ys, k = -k)
  ts2 <- cbind(ys, ys.lag)

  par(mfrow = c(3, 1), mar = c(2, 4, 1, 1))
  matplot(ts2, type = 'o', pch = 1, col = c(1, 4), 
          xlab = '', ylab = '', xlim = c(0, 120))
  abline(v = c(0, k), col = 2)
  arrows( 0, 0, k, 0, length = 0.1, col = 2)
  mtext(paste(k), side = 1, at = k, line = 1, col = 2)
  legend('topright', col = c(1, 4), bg = 'white', lty = 1,
        legend = c('原系列（周期10の正弦波＋ノイズ）', paste0(k, '次ラグ系列')))
  grid()
   acf(ys, lag.max = k, xlim = c(0, 120), ylim = c(-1, 1), main = '')
  pacf(ys, lag.max = k, xlim = c(0, 120), ylim = c(-1, 1), main = '')

  dev.off()
}
#system('bash -c 'convert -delay 100 -loop 0 ./acf_lag/acf_lag{1..25}.png ./acf_lag/acf_lag.gif'')

a <- acf(dy, lag.max = 50) # MA(q)
grid()
abline(v = fq, col = 3, lty = 3)

q <- 12 

pa <- pacf(dy, lag.max = 50) # AR(p)【注意】acfと違いグラフが1次ラグから始まる。
grid()
abline(v = fq, col = 3, lty = 3)

課題演習

n <- 24*7*2            # 時間数
t <- 1:n               # 時間
n.tr <- floor(n * 0.8) # 訓練期間数
n.te <- n - n.tr       # テスト期間数
ii.tr <- 1:n.tr        # 訓練期間インデックス
ii.te <- (n.tr+1):n    # テスト期間インデックス

set.seed(5)

y <- 100 + 0.1 * t + 0.02 * t^2 + 100 * sin(2 * pi * t / 24) + rnorm(n, sd = 5)

matplot(y, type = 'o', pch = 1, col = 2)