title: “ビジネス アナリティクス” subtitle: “演習課題” date: “2023-12-05” output: html_document: toc: yes toc_float: yes number_sections: no — — * このページは公開されるので個人情報などは記載しないこと。 * 課題が完成したらknitしPublishする。 —

# データ生成
set.seed(5)
t <- 1:(24 * 7 * 2)
y <- 100 + 0.1 * t + 0.02 * t^2 + 100 * sin(2 * pi * t / 24) + rnorm(length(t), mean = 0, sd = 5)

# データの可視化
matplot(t, y, type = 'o', pch = 1, col = 2)

fq <- 24

library(forecast)

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

# 訓練データとテストデータの分割
n.tr <- floor(length(t) * 0.8)
n.te <- length(t) - n.tr
ii.tr <- 1:n.tr
ii.te <- (n.tr + 1):length(t)

d <- 2

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

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

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

q <- 2

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

p <- 4

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

## Warning in guerrero(x, lower, upper): Guerrero's method for selecting a Box-Cox
## parameter (lambda) is given for strictly positive data.

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

## Series: ts(y[ii.tr], frequency = fq) 
## ARIMA(4,2,2)(1,1,0)[24] 
## Box Cox transformation: lambda= 0.9290708 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ma1      ma2     sar1
##       -1.2999  -0.8384  -0.4290  -0.0556  -0.3898  -0.6101  -0.5320
## s.e.   0.2203   0.1854   0.1462   0.0887   0.2095   0.2092   0.0594
## 
## sigma^2 = 21.13:  log likelihood = -719.19
## AIC=1454.39   AICc=1455.01   BIC=1482.3

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,2,2)(1,1,0)[24]
## Q* = 76.916, df = 41, p-value = 0.0005755
## 
## Model df: 7.   Total lags used: 48

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

##        MBE      MAE      MAPE     RMSE
## 1 2.120491 5.982182 0.3019825 7.755741