• このページは公開されるので個人情報などは記載しないこと。
  • 課題が完成したらknitしPublishする。

課題名 時系列分析（ARIMA）

n <- 2472 # 時間数 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)

定常化（トレンド除去）

d <- 2

dy <- diff(y[ii.tr], diff = d)

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

ACF・PACF

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

ACF

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

PACF

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

フィッティング

fit.arima <- Arima(ts(y[ii.tr], frequency = fq), lambda = ‘auto’, order = c(p, d, q), seasonal = c(1, 1, 0)) # SARIMA

c.f.

#fit.arima <- auto.arima(ts(y[ii.tr], frequency = fq), lambda = ‘auto’) fit.arima

残差分析

checkresiduals(fit.arima) # 残差分析グラフ

予測

pred <- forecast(fit.arima, h = n.te, level = 95) # グラフ matplot(t, y, type = ‘o’, pch = 16, col = 2, main = pred\(method) grid() abline(v = n.tr + 1, lty = 3, col = 3) matlines(t[ii.tr], fit.arima\)fitted, type = ‘o’, pch = 16, lty = 2, col = 4) matlines(t[ii.te], pred\(mean, type = 'o', pch = 1, lty = 2, col = 4) matlines(t[ii.te], pred\)lower, lty = 3, col = ‘gray’) matlines(t[ii.te], pred$upper, lty = 3, col = ‘gray’) legend(‘topleft’, bg = ‘white’, lty = c(1, 2, 2, 3), col = c(2, 4, 4, ‘gray’), pch = c(16, 16, 1, NA), legend = c(‘原系列’, ‘1期先予測値’, ‘予測値’, ‘95%予測区間’))

精度評価

get.accuracy <- function(yhat, y) { data.frame(MBE = mean(yhat - y), MAE = mean(abs(yhat - y)), MAPE = mean(abs((yhat - y) / y)) * 100, RMSE = sqrt(mean((yhat - y)^2)) ) }

get.accuracy(yhat = pred$mean, y = y[ii.te])